PetscSolverFeti.cc 53.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
//
// 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.


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

namespace AMDiS {

  using namespace std;

23
24
25
26
27
28
29
  // 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);
30
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
31
32

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
33
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
34
35
36
37
38
39
40
41
42
43
44
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(*(data->mat_primal_primal), x, y);
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
45
46
    //    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)
47
48
49

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

Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
55

Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
57
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
61

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

    return 0;
  }


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

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


80
    // === Restriction of the B nodes to the boundary nodes. ===
81

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

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

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
93
94

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


98
    // === K_DD - K_DI inv(K_II) K_ID ===
99

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

102
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
103
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
    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);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];

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


140
    // === Restriction of the B nodes to the boundary nodes. ===
141

142
143
144
145
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
146

147
    PetscScalar *local_b, *local_duals;
148
    VecGetArray(data->tmp_vec_b, &local_b);
149
    VecGetArray(data->tmp_vec_duals0, &local_duals);
150

151
152
    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
153
154

    VecRestoreArray(data->tmp_vec_b, &local_b);
155
156
157
158
159
160
161
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

162

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

165
166
167
168
169
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];
170

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

174
175

    // Multiply with scaled Lagrange constraint matrix.
176
177
178
179
180
181
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


182
183
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
184
185
      nComponents(-1),
      schurPrimalSolver(0)
186
187
188
189
190
191
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
192
      MSG("Create FETI-DP solver with no preconditioner!\n");
193
194
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
195
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
196
197
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
198
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
199
200
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
203
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
204
205
206
207
208

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


212
  void PetscSolverFeti::updateDofData()
213
214
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
215
216
217

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
218
   
219
220
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
221

222
223
224
225
226
      createPrimals(feSpace);      
      createDuals(feSpace);
      createLagrange(feSpace);      
      createIndexB(feSpace);
    }
227
228
229
  }


230
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
231
  {
232
    FUNCNAME("PetscSolverFeti::createPrimals()");  
233

234
235
236
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

237
238
239
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
240
    DofContainerSet& vertices = 
241
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
242
243
244
245
    TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n");
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
246

247
248
249
250

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

251
    primalDofMap.addFeSpace(feSpace);
252
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
253
254
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
255

256
    primalDofMap[feSpace].update();
257

258
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
259
260
	primalDofMap[feSpace].nRankDofs, 
	primalDofMap[feSpace].nOverallDofs);
261

262
263
264
265

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===
266
267
    
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
268

269
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
270
271
272
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
273
274
275
276
	if (primalDofMap[feSpace].isSet(it.getDofIndex())) {
	  DegreeOfFreedom d = primalDofMap[feSpace][it.getDofIndex()];
	  stdMpi.getSendData(it.getRank()).push_back(d);
	}
277

278
279
    stdMpi.updateSendDataSize();

280
281
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
282
      bool recvFromRank = false;
283
284
285
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
286
287
288
	  recvFromRank = true;
	  break;
	}
289
      }
290
291

      if (recvFromRank) 
292
	stdMpi.recv(it.getRank());
293
    }
294

295
296
    stdMpi.startCommunication();

297
298
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
299
      int i = 0;
300
301
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
302
303
304
305
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[i++];
	  primalDofMap[feSpace].insert(it.getDofIndex(), d);
	}
306
307
308
      }
    }

309
    TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size())
310
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
311
       primals.size(), primalDofMap[feSpace].size());
312
313
314
  }


315
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
316
317
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
318

319
320
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
321
322
323

    boundaryDofRanks.clear();

324
325
326
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
327
	// If DOF is not primal, i.e., its a dual node
328
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
329
330
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
331
332
	}
      }
333
	
334
335
336
337

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

338
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
339
340
341
342

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
343
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
344
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
345
346
347

    stdMpi.updateSendDataSize();

348
349
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
350
      bool recvFromRank = false;
351
      for (; !it.endDofIter(); it.nextDof()) {
352
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
353
354
355
	  recvFromRank = true;
	  break;
	}
356
      }
357
358

      if (recvFromRank)
359
	stdMpi.recv(it.getRank());
360
    }
361

362
363
    stdMpi.startCommunication();

364
365
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
366
      int i = 0;
367
      for (; !it.endDofIter(); it.nextDof())
368
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
369
370
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
371
372
373
374
    }

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

375
    dualDofMap.addFeSpace(feSpace);
376
377

    DofContainer allBoundaryDofs;
378
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
379
380

    for (DofContainer::iterator it = allBoundaryDofs.begin();
381
382
383
	 it != allBoundaryDofs.end(); ++it)
      if (primalDofMap[feSpace].isSet(**it) == false)
	dualDofMap[feSpace].insertRankDof(**it);
384

385
    dualDofMap[feSpace].update(false);
386

387
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
388
389
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
390
391
392
  }

  
393
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
394
395
396
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

397
398
399
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

400
    lagrangeMap.addFeSpace(feSpace);
401
402
403
404
405

    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
406
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
407
	int degree = boundaryDofRanks[it->first].size();
408
409
410
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
411
412
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
    lagrangeMap[feSpace].update();
413
    
414
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
415
416
	lagrangeMap[feSpace].nRankDofs, 
	lagrangeMap[feSpace].nOverallDofs);
417
418


419
    // === Communicate lagrangeMap to all other ranks. ===
420
421

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
422

423
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
424
	 !it.end(); it.nextRank()) {
425
      for (; !it.endDofIter(); it.nextDof())
426
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
427
	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
428
	  stdMpi.getSendData(it.getRank()).push_back(d);
429
	}
430
    }
431

432
433
    stdMpi.updateSendDataSize();

434
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
435
	 !it.end(); it.nextRank()) {
436
      bool recvData = false;
437
      for (; !it.endDofIter(); it.nextDof())
438
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
439
440
441
442
443
	  recvData = true;
	  break;
	}
	  
      if (recvData)
444
	stdMpi.recv(it.getRank());
445
446
447
448
    }

    stdMpi.startCommunication();

449
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
450
	 !it.end(); it.nextRank()) {
451
      int counter = 0;
452
453
454
      for (; !it.endDofIter(); it.nextDof()) {
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
455
	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
456
457
	}
      }
458
    }
459
460
461
  }


462
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
463
  {
464
    FUNCNAME("PetscSolverFeti::createIndexB()");
465

466
    localDofMap.addFeSpace(feSpace);
467
    DOFAdmin* admin = feSpace->getAdmin();
468
469
470
471

    // === 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.                                 ===
472
473
474

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
475
476
      if (admin->isDofFree(i) == false && 
	  primalDofMap[feSpace].isSet(i) == false && 
477
478
479
480
	  dualDofMap[feSpace].isSet(i) == false) {
	localDofMap[feSpace].insertRankDof(i);
	nLocalInterior++;
      }
481
482
    }

483
484
    localDofMap[feSpace].update(false);

485
486
    TEST_EXIT_DBG(nLocalInterior + 
		  primalDofMap[feSpace].size() + 
487
		  dualDofMap[feSpace].size() == 
488
489
490
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

491
492
    // === And finally, add the global indicies of all dual nodes. ===

493
494
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
495
496
497
498
499
500
      localDofMap[feSpace].insertRankDof(it->first);

    localDofMap[feSpace].update(false);

    dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + 
				  nLocalInterior);
501
502
503
  }


504
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
505
506
507
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

508
509
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

510
511
    // === Create distributed matrix for Lagrange constraints. ===

512
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
513
		    lagrangeMap.getRankDofs(feSpaces),
514
515
516
		    localDofMap.getRankDofs(feSpaces),
		    lagrangeMap.getOverallDofs(feSpaces),
		    localDofMap.getOverallDofs(feSpaces),	
517
518
519
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

520
521
522
523
524
525
526
    // === 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.                                                     ===

527
528
529
530
    DofMapping &dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");
531

532
      // Global index of the first Lagrange constriant for this node.
533
      int index = lagrangeMap[feSpace][it->first];
534
      // Copy set of all ranks that contain this dual node.
535
536
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
537
      // Number of ranks that contain this dual node.
538
539
540
541
      int degree = W.size();

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
542
543
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
544
545
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
546
	      int colIndex = it->second * nComponents + k;
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
	    }
	  }

	  index++;
	}
      }
    }

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


563
  void PetscSolverFeti::createSchurPrimalKsp()
564
565
566
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

567
568
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
574
575
576
577
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
      VecCreateMPI(PETSC_COMM_WORLD, 
578
		   localDofMap[feSpace].nRankDofs * nComponents, localDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
579
580
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
581
		   primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
582
583
584
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
585
586
		     primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
		     primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
587
588
589
590
591
592
593
594
595
596
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
    } else {
597
598
      MSG("Create direct schur primal solver!\n");

599
600
601
602
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

603
604
      int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
      int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
605
606
      int nRowsRankB = localDofMap[feSpace].nRankDofs * nComponents;
      int nRowsOverallB = localDofMap[feSpace].nOverallDofs * nComponents;
607

Thomas Witkowski's avatar
Thomas Witkowski committed
608
609
610
611
612
613
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
614

615
616
617
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
618

619
620
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

621
622
      for (int i = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++) {
	PetscInt row = localDofMap[feSpace].rStartDofs * nComponents + i;
623
624
625
626
627
628
	MatGetRow(mat_b_primal, row, &nCols, &cols, &values);

	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
629
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
630
631
632
633
634
635
      }

      int maxLocalPrimal = mat_b_primal_cols.size();
      mpi::globalMax(maxLocalPrimal);

      TEST_EXIT(mat_b_primal_cols.size() ==
636
		primalDofMap[feSpace].size() * nComponents)("Should not happen!\n");
637
638
639
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
640
	VecCreateSeq(PETSC_COMM_SELF, localDofMap[feSpace].nRankDofs * nComponents, &tmpVec);
641
642
643
644
645
646
647
648
649
650

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

       	KSPSolve(ksp_b, tmpVec, tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
651
	PetscScalar *tmpValues;
652
	VecGetArray(tmpVec, &tmpValues);
653
	for (int i  = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
654
	  MatSetValue(matBPi, 
655
		      localDofMap[feSpace].rStartDofs * nComponents + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
658
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
659
660
661
662
663
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
667
668
669
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
670
      MatDestroy(&matBPi);
671
672
673
674
675
676
677
678
679
680

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

      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_primal_primal, 
		      mat_primal_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
681

682
683
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
684
    }
685
686
687
688
689
690
691
  }


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

692
693
694
695
696
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
697

698
699
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
700

701
702
703
704
705
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
706
707
708
709
710
711
712
  }


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

713
714
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

715
716
    // === Create FETI-DP solver object. ===

717
718
719
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
720
    fetiData.fetiSolver = this;
721
    fetiData.ksp_schur_primal = &ksp_schur_primal;
722

723
    VecCreateMPI(PETSC_COMM_WORLD, 
724
		 localDofMap[feSpace].nRankDofs * nComponents, localDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
727
728
		 lagrangeMap[feSpace].nRankDofs * nComponents, 
		 lagrangeMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
729
		 &(fetiData.tmp_vec_lagrange));
730
    VecCreateMPI(PETSC_COMM_WORLD, 
731
		 primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
732
		 &(fetiData.tmp_vec_primal));
733
734

    MatCreateShell(PETSC_COMM_WORLD,
735
736
737
738
		   lagrangeMap[feSpace].nRankDofs * nComponents, 
		   lagrangeMap[feSpace].nRankDofs * nComponents,
		   lagrangeMap[feSpace].nOverallDofs * nComponents, 
		   lagrangeMap[feSpace].nOverallDofs * nComponents,
739
		   &fetiData, &mat_feti);
740
741
742
743
744
745
746
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


    KSPCreate(PETSC_COMM_WORLD, &ksp_feti);
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
    KSPSetFromOptions(ksp_feti);
747
748


749
    // === Create FETI-DP preconditioner object. ===
750

751
752
753
754
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
755

756
757
758
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
761
762
763
764
765
766
767
768
769
770
      KSPSetOptionsPrefix(ksp_interior, "solver_interior_");
      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;
      
771
      VecCreateMPI(PETSC_COMM_WORLD, 
772
		   localDofMap[feSpace].nRankDofs * nComponents, localDofMap[feSpace].nOverallDofs * nComponents,
773
		   &(fetiDirichletPreconData.tmp_vec_b));      
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
      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));
      
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

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

789
      VecCreateMPI(PETSC_COMM_WORLD, 
790
		   localDofMap[feSpace].nRankDofs * nComponents, localDofMap[feSpace].nOverallDofs * nComponents,
791
		   &(fetiLumpedPreconData.tmp_vec_b));
792
793
794
795
796
797
798
799
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
800
801
      break;
    default:
802
803
      break;
    }
804
805
806
807
808
809
810
  }
  

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

811
812
    // === Destroy FETI-DP solver object. ===

813
814
815
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
816
    fetiData.fetiSolver = NULL;
817
    fetiData.ksp_schur_primal = PETSC_NULL;
818

Thomas Witkowski's avatar
Thomas Witkowski committed
819
820
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
821
    VecDestroy(&fetiData.tmp_vec_primal);
822
823
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
824
825


826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
    // === Destroy FETI-DP preconditioner object. ===

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPDestroy(&ksp_interior);

      fetiDirichletPreconData.mat_lagrange_scaled = NULL;
      fetiDirichletPreconData.mat_interior_interior = NULL;
      fetiDirichletPreconData.mat_duals_duals = NULL;
      fetiDirichletPreconData.mat_interior_duals = NULL;
      fetiDirichletPreconData.mat_duals_interior = NULL;
      fetiDirichletPreconData.ksp_interior = NULL;
      
      VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
      MatDestroy(&mat_lagrange_scaled);
      break;

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

      VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
      break;
854
855
    default:
      break;
856
    }
857
858
859
860
861
862
863
864
865
  }


  void PetscSolverFeti::recoverSolution(Vec &vec_sol_b,
					Vec &vec_sol_primal,
					SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::recoverSolution()");

866
867
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

868
    // === Get local part of the solution for B variables. ===
869
870
871
872
873

    PetscScalar *localSolB;
    VecGetArray(vec_sol_b, &localSolB);


874
875
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
876
877
    
    vector<PetscInt> globalIsIndex, localIsIndex;
878
879
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
880
881
882

    {
      int counter = 0;
883
884
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
	for (int i = 0; i < nComponents; i++) {
	  globalIsIndex.push_back(it->second * nComponents + i);
	  localIsIndex.push_back(counter++);
	}
      }
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);

    ISCreateGeneral(PETSC_COMM_SELF, 
		    localIsIndex.size(), 
		    &(localIsIndex[0]),
		    PETSC_USE_POINTER,
		    &localIs);

    Vec local_sol_primal;
    VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_primal);

    VecScatter primalScatter;
    VecScatterCreate(vec_sol_primal, globalIs, local_sol_primal, localIs, &primalScatter);
    VecScatterBegin(primalScatter, vec_sol_primal, local_sol_primal, 
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(primalScatter, vec_sol_primal, local_sol_primal, 
		  INSERT_VALUES, SCATTER_FORWARD);

915
916
917
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
918
919
920
921

    PetscScalar *localSolPrimal;
    VecGetArray(local_sol_primal, &localSolPrimal);

922
    // === And copy from PETSc local vectors to the DOF vectors. ===
923
924
925
926

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

927
928
      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
	   it != localDofMap[feSpace].getMap().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
929
	int petscIndex = it->second * nComponents + i;
930
931
932
933
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
934
935
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
936
937
938
939
940
941
942
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
943
    VecDestroy(&local_sol_primal);
944
945
946
  }


947
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
948
949
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
950

951
952
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

953
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);