PetscSolver.cc 4.1 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.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include "parallel/PetscSolver.h"
14

15
16
17

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
18
19
  using namespace std;

20
21
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
  {    
22
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
23
      cout << "[0]  Petsc-Iteration " << iter << ": " << rnorm << endl;
24
25
26
27

    return 0;
  }

28

29
  PetscSolver* PetscSolver::getPetscSolver(Mat& petscMatrix, MeshDistributor *dist, int n)
30
  {
31
    FUNCNAME("PetscSolver::getPetscSolver()");
32

33
34
    string name("");
    GET_PARAMETER(0, "parallel->solver", &name);
35
    
36
37
38
39
40
41
    if (name == "petsc-schur") {
      return new PetscSolverSchur(petscMatrix, dist, n);
    } else if (name == "petsc" || name == "") {
      return new PetscSolverGlobalMatrix(petscMatrix, dist, n);
    } else {
      ERROR_EXIT("No parallel solver %s available!\n", name.c_str());
42
43
44
45
    }
  }


46
  void PetscSolverGlobalMatrix::providePetscSolver(KSP &solver, PC &pc)
47
  {
48
49
50
51
52
53
54
55
56
57
58
    FUNCNAME("PetscSolverGlobalMatrix::providePetscProblemStat()");
    
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPGetPC(solver, &pc);
    
    KSPSetOperators(solver, matrix, matrix, SAME_NONZERO_PATTERN); 
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
    KSPSetFromOptions(solver);
    PCSetFromOptions(pc);
59
60
61
  }


62
  void PetscSolverSchur::providePetscSolver(KSP &solver, PC &pc)
63
  {
64
65
    FUNCNAME("PetscSolverSchur::providePetscProblemStat()");
    
Thomas Witkowski's avatar
Thomas Witkowski committed
66
    typedef map<int, DofContainer> RankToDofContainer;
67
68
69
70
71
72
73
74
75
76
    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);
77
	  
78
79
	for (int i = 0; i < nComponents; i++)
	  boundaryDofsSet.insert(meshDistributor->mapLocalToGlobal(**dofIt) * nComponents + i);
80
      }
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
      
    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());
105
106
107

    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPGetPC(solver, &pc);
108
109
110
111
112
113
114
115
116
117
118
119
120
121
    PCSetType(pc, PCFIELDSPLIT);
      
    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);
     
     
    KSPSetOperators(solver, matrix, matrix, SAME_NONZERO_PATTERN); 
122
123
124
125
126
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetFromOptions(solver);
    PCSetFromOptions(pc);
  }

127
}