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


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

namespace AMDiS {

  using namespace std;

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


31 32 33 34 35 36 37
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

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

40 41 42 43
    MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
44 45 46 47 48 49 50 51 52
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


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

55 56
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
57

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

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

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

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

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

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

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

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

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

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

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

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

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

91 92 93 94
    return 0;
  }


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

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

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


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

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

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

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

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


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

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

132
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
133
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
134 135 136 137 138 139 140 141 142 143
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

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


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

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

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

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


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

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

157 158 159 160 161 162 163 164 165 166 167 168 169 170
    return 0;
  }


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

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


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

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

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

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

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


    // === K_DD ===

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

196

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

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

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

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

209 210

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

    return 0;
  }


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

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

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

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

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


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

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

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

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

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

282 283 284 285
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
286

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

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

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


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

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

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

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

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

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

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

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

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

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

342 343
      createPrimals(feSpace);  

344 345 346
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
347

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

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

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

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

365 366 367
    lagrangeMap.update();


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

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

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

380
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
381 382 383 384 385 386 387 388 389
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

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

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

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

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

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

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


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

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

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

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
448 449 450 451 452 453 454
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
455
    }
456

457 458 459 460

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

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


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

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

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

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

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

  
493 494 495 496 497 498 499 500 501 502 503 504
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (feSpace != fullInterface)
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
505 506
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
507
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	interfaceDofMap[feSpace].insertNonRankDof(**it);
509 510 511
  }


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

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

519 520
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
521 522 523
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
524
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
525

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

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
537
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
538
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
539
	  else
540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
555 556 557 558 559 560
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
561

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
562 563 564 565 566 567 568 569 570 571 572 573 574 575
    if (dualDofMap[feSpace].nLocalDofs == 0)
      return;


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

    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(feSpace, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

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

583 584 585 586

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

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

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

    stdMpi.updateSendDataSize();

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

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

617 618
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
633

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

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


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

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

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
661

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

671 672 673 674 675 676 677 678 679
	  if (fetiPreconditioner != FETI_NONE)
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
680 681 682

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

689 690 691 692 693 694 695 696 697 698
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
699 700 701
  }


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

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

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

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

717 718 719 720 721 722 723
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
724
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
725
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
726

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
741 742 743 744
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
745 746
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
747

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

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

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


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

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

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

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

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

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

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

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

820

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

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

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

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

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