PetscSolverSchur.cc 2.82 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
//
// 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 <vector>
#include <set>

#include "parallel/PetscSolverSchur.h"

namespace AMDiS {

  void PetscSolverSchur::providePetscSolver(KSP &solver, PC &pc)
  {
    FUNCNAME("PetscSolver::providePetscSolver()");

    typedef map<int, DofContainer> RankToDofContainer;
    typedef map<DegreeOfFreedom, bool> DofIndexToBool;

    std::set<DegreeOfFreedom> boundaryDofsSet;
    std::set<DegreeOfFreedom> boundaryLocalDofs;
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator rankIt = sendDofs.begin(); 
	 rankIt != sendDofs.end(); ++rankIt)
      for (DofContainer::iterator dofIt = rankIt->second.begin(); 
	   dofIt != rankIt->second.end(); ++dofIt) {
	boundaryLocalDofs.insert(**dofIt);

	for (int i = 0; i < nComponents; i++)
	  boundaryDofsSet.insert(meshDistributor->mapLocalToGlobal(**dofIt) * nComponents + i);
      }

    vector<DegreeOfFreedom> boundaryDofs(boundaryDofsSet.begin(), 
					 boundaryDofsSet.end());

    std::set<DegreeOfFreedom> otherBoundaryLocalDofs;
    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator rankIt = recvDofs.begin(); 
	 rankIt != recvDofs.end(); ++rankIt)
      for (DofContainer::iterator dofIt = rankIt->second.begin(); 
	   dofIt != rankIt->second.end(); ++dofIt)
	otherBoundaryLocalDofs.insert(**dofIt);

    std::set<DegreeOfFreedom> interiorDofsSet;
    DofIndexToBool& isRankDof = meshDistributor->getIsRankDof();
    for (DofIndexToBool::iterator dofIt = isRankDof.begin(); 
	 dofIt != isRankDof.end(); ++dofIt)
      if (dofIt->second && 
	  boundaryLocalDofs.count(dofIt->first) == 0 && 
	  otherBoundaryLocalDofs.count(dofIt->first) == 0)
	for (int i = 0; i < nComponents; i++)	    
	  interiorDofsSet.insert(meshDistributor->mapLocalToGlobal(dofIt->first) * nComponents + i);

    vector<DegreeOfFreedom> interiorDofs(interiorDofsSet.begin(), 
					 interiorDofsSet.end());

    IS interiorIs;
    ISCreateGeneral(PETSC_COMM_WORLD, interiorDofs.size(), &(interiorDofs[0]), 
		    PETSC_COPY_VALUES, &interiorIs);
    PCFieldSplitSetIS(pc, "interior", interiorIs);

    IS boundaryIs;
    ISCreateGeneral(PETSC_COMM_WORLD, boundaryDofs.size(), &(boundaryDofs[0]), 
		    PETSC_COPY_VALUES, &boundaryIs);
    PCFieldSplitSetIS(pc, "boundary", boundaryIs);

    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPGetPC(solver, &pc);
    PCSetType(pc, PCFIELDSPLIT);

    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetFromOptions(solver);
    PCSetFromOptions(pc);
  }

}