Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

PetscSolverFeti.cc 40.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 "parallel/SubDomainSolver.h"
18
#include "io/VtkWriter.h"
19 20 21 22 23

namespace AMDiS {

  using namespace std;

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

33 34 35 36
    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);
37 38 39 40 41 42 43 44 45
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

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

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

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

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

    return 0;
  }


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

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


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

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

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

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

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


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

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

104
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
105
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
106 107 108 109 110 111 112 113 114 115
    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);

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

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


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

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

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

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

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


    // === K_DD ===

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

166

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

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

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

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

179 180

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

    return 0;
  }


187 188
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
189
      schurPrimalSolver(0),
190 191
      multiLevelTest(false),
      subDomainSolver(NULL)
192 193 194 195 196 197
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
217 218 219
  }


220
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
221
  {
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
    FUNCNAME("PetscSolverFeti::initialize()");

    if (subDomainSolver == NULL)
      subDomainSolver = new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);

    primalDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), 
		      true, true);
    dualDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		    false, false);
    lagrangeMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		     true, true);
    localDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(),
		     false, false);
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.init(mpiComm, feSpaces,  meshDistributor->getFeSpaces(),
			  false, false);
  }


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

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

249 250 251 252 253 254
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
    interiorDofMap.clear();

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

258 259
      createPrimals(feSpace);      
      createDuals(feSpace);
260
      createLagrange(feSpace);      
261 262
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
263

264 265 266 267 268
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
269 270 271 272
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
273
    if (fetiPreconditioner != FETI_NONE)
274
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
275 276 277

    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
278
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
279

280
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
281 282
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
283
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
284
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
285

286
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
287
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
288 289 290 291

      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
292
    }
293 294 295


    // If multi level test, inform sub domain solver about coarse space.
296
    subDomainSolver->setDofMapping(&primalDofMap, &localDofMap);
297 298 299
  }


300
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
301
  {
302
    FUNCNAME("PetscSolverFeti::createPrimals()");  
303

304 305 306
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

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

317 318 319 320

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

321
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
322 323
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
324 325
      else
  	primalDofMap[feSpace].insert(*it);
326 327 328
  }


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

333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);

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

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

349 350
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
351

352
    boundaryDofRanks[feSpace].clear();
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 359
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
360 361
	}
      }
362
	
363 364 365 366

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

367
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
368 369 370 371

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
372
	if (!isPrimal(feSpace, it.getDofIndex()))
373
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
374 375 376

    stdMpi.updateSendDataSize();

377 378
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
379
      bool recvFromRank = false;
380
      for (; !it.endDofIter(); it.nextDof()) {
381
	if (!isPrimal(feSpace, it.getDofIndex())) {
382 383 384
	  recvFromRank = true;
	  break;
	}
385
      }
386 387

      if (recvFromRank)
388
	stdMpi.recv(it.getRank());
389
    }
390

391 392
    stdMpi.startCommunication();

393 394
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
395
      int i = 0;
396
      for (; !it.endDofIter(); it.nextDof())
397
	if (!isPrimal(feSpace, it.getDofIndex()))
398
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
399
	    stdMpi.getRecvData(it.getRank())[i++];
400 401 402
    }


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

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


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

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

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

431
    int nLocalInterior = 0;    
432
    for (int i = 0; i < admin->getUsedSize(); i++) {
433
      if (admin->isDofFree(i) == false && 
434
	  isPrimal(feSpace, i) == false &&
435
	  dualDofMap[feSpace].isSet(i) == false) {
436 437
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

438
	if (fetiPreconditioner != FETI_NONE)
439 440 441
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
442
      }
443
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
444
    
445 446
    // === And finally, add the global indicies of all dual nodes. ===

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


453
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
454 455 456
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

457 458
    // === Create distributed matrix for Lagrange constraints. ===

459
    MatCreateMPIAIJ(mpiComm,
460 461
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
462 463 464
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

465 466 467 468 469 470 471
    // === 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
472 473 474
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
      map<DegreeOfFreedom, MultiIndex> &dualMap = 
	dualDofMap[feSpaces[k]].getMap();
475

Thomas Witkowski's avatar
Thomas Witkowski committed
476 477
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	   it != dualMap.end(); ++it) {
478
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
479 480 481 482 483 484
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);
	
	// Copy set of all ranks that contain this dual node.
485 486
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
487 488 489 490 491 492
	// 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) {
493
	      int colIndex = localDofMap.getMatIndex(k, it->first);
494
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
495
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
496
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
497
	    index++;	      
498 499 500 501 502 503 504 505 506 507
	  }
	}
      }
    }

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


508
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
509 510 511
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

515
      schurPrimalData.subSolver = subDomainSolver;
Thomas Witkowski's avatar
Thomas Witkowski committed
516
      
517
      VecCreateMPI(mpiComm, 
518
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
519
		   &(schurPrimalData.tmp_vec_b));
520
      VecCreateMPI(mpiComm, 
521
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
522 523
		   &(schurPrimalData.tmp_vec_primal));

524
      MatCreateShell(mpiComm,
525 526
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
527 528 529 530 531
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
532
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
533
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
534 535
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
536 537
      KSPSetFromOptions(ksp_schur_primal);
    } else {
538 539
      MSG("Create direct schur primal solver!\n");

540 541
      double wtime = MPI::Wtime();

542 543
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
544 545
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
546

Thomas Witkowski's avatar
Thomas Witkowski committed
547
      Mat matBPi;
548
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
549 550 551 552
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
553

554 555 556
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
557

558 559
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

560 561
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
562
	MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
563 564 565 566 567

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

568
	MatRestoreRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
569 570
      }

571 572
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
573
	("Should not happen!\n");
574

575 576 577
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
578
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
579 580 581 582 583 584 585 586

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

587
	subDomainSolver->solve(tmpVec, tmpVec);
588

Thomas Witkowski's avatar
Thomas Witkowski committed
589
	PetscScalar *tmpValues;
590
	VecGetArray(tmpVec, &tmpValues);
591
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
592
	  MatSetValue(matBPi, 
593
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
594 595 596
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
597 598 599 600 601
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
602 603
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
604 605

      MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
606
      MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
607
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
608 609

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
610
      MatDestroy(&matBPi);
611 612

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

616
      KSPCreate(mpiComm, &ksp_schur_primal);
617
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
618
		      SAME_NONZERO_PATTERN);
619 620 621 622 623 624
      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);
625
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
626

627 628
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
629
    }
630 631 632 633 634 635 636
  }


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

637
    if (schurPrimalSolver == 0) {
638
      schurPrimalData.subSolver = NULL;
639

640 641 642
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
643 644 645
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
646 647 648
  }


649
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
650 651 652
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

653 654
    // === Create FETI-DP solver object. ===

655
    fetiData.mat_lagrange = &mat_lagrange;
656
    fetiData.subSolver = subDomainSolver;
657
    fetiData.ksp_schur_primal = &ksp_schur_primal;
658

659
    VecCreateMPI(mpiComm, 
660
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
661
		 &(fetiData.tmp_vec_b));
662
    VecCreateMPI(mpiComm,
663
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
664
		 &(fetiData.tmp_vec_lagrange));
665
    VecCreateMPI(mpiComm, 
666
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
667
		 &(fetiData.tmp_vec_primal));
668

669
    MatCreateShell(mpiComm,
670 671
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
672
		   &fetiData, &mat_feti);
673 674 675
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


676
    KSPCreate(mpiComm, &ksp_feti);
677
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
678 679 680
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
681
    KSPSetFromOptions(ksp_feti);
682 683


684
    // === Create FETI-DP preconditioner object. ===
685

686 687 688 689
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
690

691 692 693
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
694 695
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
696 697 698 699 700 701
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
702 703 704 705 706 707 708 709 710
      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;
      
711
      VecCreateMPI(mpiComm, 
712
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
713
		   &(fetiDirichletPreconData.tmp_vec_b));      
714 715 716 717 718 719
      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));
720 721 722 723 724 725 726 727 728 729 730 731 732

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

733 734 735 736 737 738 739 740 741 742 743
      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;

Thomas Witkowski's avatar
Thomas Witkowski committed
744 745 746 747 748 749 750 751 752 753 754 755
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

756
      VecCreateMPI(mpiComm, 
757 758
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
759
		   &(fetiLumpedPreconData.tmp_vec_b));
760 761 762 763
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
764 765 766 767 768 769

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
770 771
      break;
    default:
772 773
      break;
    }
774 775 776 777 778 779 780
  }
  

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

781 782
    // === Destroy FETI-DP solver object. ===

783
    fetiData.mat_lagrange = PETSC_NULL;
784
    fetiData.subSolver = NULL;
785
    fetiData.ksp_schur_primal = PETSC_NULL;
786

Thomas Witkowski's avatar
Thomas Witkowski committed
787 788
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
789
    VecDestroy(&fetiData.tmp_vec_primal);
790 791
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
792 793


794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821
    // === 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;
822 823
    default:
      break;
824
    }
825 826 827 828 829 830 831 832 833
  }


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

834
    // === Get local part of the solution for B variables. ===
835 836 837 838 839

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


840 841
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
842
    
Thomas Witkowski's avatar
Thomas Witkowski committed
843 844
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);

845
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
846 847
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
848 849

    {