PetscSolverFeti.cc 45 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
      createPrimals(feSpace);      
      createDuals(feSpace);
224
      createLagrange(feSpace);      
225 226
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291

    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();

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

      MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
      MSG("nRankDuals = %d   nOverallDuals = %d\n",
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);

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

    // === Communicate lagrangeMap to all other ranks. ===

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

      StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
      
      for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	   !it.end(); it.nextRank()) {
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()].local;
	    stdMpi.getSendData(it.getRank()).push_back(d);
	  }
      }
      
      stdMpi.updateSendDataSize();
      
      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	   !it.end(); it.nextRank()) {
	bool recvData = false;
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    recvData = true;
	    break;
	  }
	
	if (recvData)
	  stdMpi.recv(it.getRank());
      }
      
      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	   !it.end(); it.nextRank()) {
	int counter = 0;
	for (; !it.endDofIter(); it.nextDof()) {
	  if (!isPrimal(feSpace, it.getDofIndex())) {
	    DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
	    lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
	  }
	}
      }
    }
292 293 294
  }


295
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
296
 {
297
    FUNCNAME("PetscSolverFeti::createPrimals()");  
298

299 300 301
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

302 303 304
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
305
    DofContainerSet& vertices = 
306
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
307 308 309 310
    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);
311

312 313 314 315

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

316
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
317 318
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
319 320
      else
  	primalDofMap[feSpace].insert(*it);
321

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
322 323 324
    primalDofMap[feSpace].setOverlap(true);
    primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
  				     meshDistributor->getRecvDofs());
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
    primalDofMap[feSpace].update();
327 328 329
  }


330
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
331 332
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
333

334 335
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
336 337 338

    boundaryDofRanks.clear();

339 340 341
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
342
	if (!isPrimal(feSpace, it.getDofIndex())) {
343 344
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
345 346
	}
      }
347
	
348 349 350 351

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

352
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
353 354 355 356

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
357
	if (!isPrimal(feSpace, it.getDofIndex()))
358
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
359 360 361

    stdMpi.updateSendDataSize();

362 363
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
364
      bool recvFromRank = false;
365
      for (; !it.endDofIter(); it.nextDof()) {
366
	if (!isPrimal(feSpace, it.getDofIndex())) {
367 368 369
	  recvFromRank = true;
	  break;
	}
370
      }
371 372

      if (recvFromRank)
373
	stdMpi.recv(it.getRank());
374
    }
375

376 377
    stdMpi.startCommunication();

378 379
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
380
      int i = 0;
381
      for (; !it.endDofIter(); it.nextDof())
382
	if (!isPrimal(feSpace, it.getDofIndex()))
383 384
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
385 386 387 388 389
    }

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

    DofContainer allBoundaryDofs;
390
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
391 392

    for (DofContainer::iterator it = allBoundaryDofs.begin();
393
	 it != allBoundaryDofs.end(); ++it)
394
      if (!isPrimal(feSpace, **it))
395
	dualDofMap[feSpace].insertRankDof(**it);
396

Thomas Witkowski's avatar
Thomas Witkowski committed
397
    dualDofMap[feSpace].update();
398 399 400
  }

  
401
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
402 403 404
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

405 406 407
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

408
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
409 410
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
411
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
412
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
413
	int degree = boundaryDofRanks[it->first].size();
414 415 416
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
417
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
418
    lagrangeMap[feSpace].update();
419 420 421
  }


422
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
423
  {
424
    FUNCNAME("PetscSolverFeti::createIndexB()");
425

426
    DOFAdmin* admin = feSpace->getAdmin();
427 428 429 430

    // === 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.                                 ===
431 432 433

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
434
      if (admin->isDofFree(i) == false && 
435
	  isPrimal(feSpace, i) == false &&
436
	  dualDofMap[feSpace].isSet(i) == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
437
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
438 439
	nLocalInterior++;
      }
440
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
441
    
Thomas Witkowski's avatar
Thomas Witkowski committed
442
    TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + 
443
		  dualDofMap[feSpace].size() == 
444 445 446
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

447 448
    // === And finally, add the global indicies of all dual nodes. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
450
	 it != dualDofMap[feSpace].getMap().end(); ++it)
451 452
      localDofMap[feSpace].insertRankDof(it->first);

Thomas Witkowski's avatar
Thomas Witkowski committed
453
    localDofMap[feSpace].update();
454 455 456
  }


457
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
458 459 460
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

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

463 464
    // === Create distributed matrix for Lagrange constraints. ===

465
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
466 467
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
468 469 470
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

471 472 473 474 475 476 477
    // === 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
478 479
    map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
480 481 482 483
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");

      // Global index of the first Lagrange constriant for this node.
Thomas Witkowski's avatar
Thomas Witkowski committed
484
      int index = lagrangeMap[feSpace][it->first].local;
485 486 487 488 489 490 491 492 493 494 495 496
      // Copy set of all ranks that contain this dual node.
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
      // Number of ranks that contain this dual node.
      int degree = W.size();

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
Thomas Witkowski's avatar
Thomas Witkowski committed
497 498 499 500
	      //	      int colIndex = it->second * nComponents + k;
	      int colIndex = 
		(localDofMap[feSpace].rStartDofs +
		 localDofMap[feSpace][it->first].local) * nComponents + k;
501
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
502 503
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
504 505
	    }
	  }
506 507

	  index++;
508 509 510 511 512 513 514 515 516
	}
      }
    }

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


517
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
518 519 520
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
526 527 528 529 530 531
      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, 
532
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
533 534
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
535
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
536 537 538
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
539 540
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
541 542 543 544 545 546 547 548 549 550
		     &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 {
551 552
      MSG("Create direct schur primal solver!\n");

553 554 555 556
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

557 558
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
559 560
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
561

Thomas Witkowski's avatar
Thomas Witkowski committed
562 563 564 565 566 567
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
568

569 570 571
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
572

573 574
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

575 576
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
577 578 579 580 581 582
	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
583
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
584 585 586 587 588 589
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
590 591
		primalDofMap[feSpace].size() * nComponents)
	("Should not happen!\n");
592 593 594
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
595
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
596 597 598 599 600 601 602 603 604 605

 	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
606
	PetscScalar *tmpValues;
607
	VecGetArray(tmpVec, &tmpValues);
608
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	  MatSetValue(matBPi, 
610
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
611 612 613
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
614 615 616 617 618
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
619 620 621
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
622 623 624
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
625
      MatDestroy(&matBPi);
626 627 628 629 630 631 632 633 634 635

      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
636

637 638
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
639
    }
640 641 642 643 644 645 646
  }


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

647 648 649 650 651
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
652

653 654
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
655

656 657 658 659 660
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
661 662 663
  }


664
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
665 666 667
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

668 669
    // === Create FETI-DP solver object. ===

670 671 672
    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
673
    fetiData.fetiSolver = this;
674
    fetiData.ksp_schur_primal = &ksp_schur_primal;
675

676
    VecCreateMPI(PETSC_COMM_WORLD, 
677
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
678 679
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
680
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
681
		 &(fetiData.tmp_vec_lagrange));
682
    VecCreateMPI(PETSC_COMM_WORLD, 
683
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
684
		 &(fetiData.tmp_vec_primal));
685 686

    MatCreateShell(PETSC_COMM_WORLD,
687 688
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
689
		   &fetiData, &mat_feti);
690 691 692 693 694 695 696
    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);
697 698


699
    // === Create FETI-DP preconditioner object. ===
700

701 702 703 704
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
705

706 707 708
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
709 710
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
711 712 713 714 715 716 717 718 719 720
      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;
      
721
      VecCreateMPI(PETSC_COMM_WORLD, 
722
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
723
		   &(fetiDirichletPreconData.tmp_vec_b));      
724 725 726 727 728 729
      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));
730 731 732 733 734 735 736 737 738 739 740 741
      
      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;

742
      VecCreateMPI(PETSC_COMM_WORLD, 
743 744
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
745
		   &(fetiLumpedPreconData.tmp_vec_b));
746 747 748 749
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
750 751 752 753 754 755

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
756 757
      break;
    default:
758 759
      break;
    }
760 761 762 763 764 765 766
  }
  

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

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

769 770 771
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
772
    fetiData.fetiSolver = NULL;
773
    fetiData.ksp_schur_primal = PETSC_NULL;
774

Thomas Witkowski's avatar
Thomas Witkowski committed
775 776
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
777
    VecDestroy(&fetiData.tmp_vec_primal);
778 779
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
780 781


782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809
    // === 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;
810 811
    default:
      break;
812
    }
813 814 815 816 817 818 819 820 821
  }


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

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

824
    // === Get local part of the solution for B variables. ===
825 826 827 828 829

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


830 831
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
832 833
    
    vector<PetscInt> globalIsIndex, localIsIndex;
834 835
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
836 837 838

    {
      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
839
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
840
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
841
	for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
842
	  globalIsIndex.push_back(it->second.local * nComponents + i);
843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
	  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);

871 872 873
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
874 875 876 877

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

878
    // === And copy from PETSc local vectors to the DOF vectors. ===
879 880 881 882

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

Thomas Witkowski's avatar
Thomas Witkowski committed
883
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpace].getMap().begin();
884
	   it != localDofMap[feSpace].getMap().end(); ++it) {
885
	int petscIndex = localDofMap.mapLocal(it->first, i);
886 887 888 889
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
890
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
891
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
892 893 894 895 896 897 898
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
899
    VecDestroy(&local_sol_primal);
900 901 902
  }


903
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
904 905
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
906

907 908
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

909
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
910
    nComponents = mat->getSize();
911 912 913

    // === Create all sets and indices. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
914 915 916 917
    primalDofMap.init(mpiComm, feSpaces, true);
    dualDofMap.init(mpiComm, feSpaces, false);
    lagrangeMap.init(mpiComm, feSpaces, true);
    localDofMap.init(mpiComm, feSpaces, false);
918

919 920
    updateDofData();

921 922
    // === Create matrices for the FETI-DP method. ===

923 924
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
925 926 927
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
    int nRowsDual = dualDofMap.getRankDofs();
928 929
    int nRowsInterior = nRowsRankB - nRowsDual;

Thomas Witkowski's avatar
Thomas Witkowski committed
930 931
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
932 933

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
934 935
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
936
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
937 938

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
939 940
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
941
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
942 943

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
944 945
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
946
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
947

948

949 950 951 952
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
953
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
954 955 956 957
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
958
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
959 960 961
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,