PetscSolverNavierStokes.cc 5.01 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
//
// 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/PetscSolverNavierStokes.h"
15
#include "parallel/PetscHelper.h"
16
17
18
19
20

namespace AMDiS {

  using namespace std;

21

22
  PetscErrorCode pcSchurShell(PC pc, Vec x, Vec y)
Thomas Witkowski's avatar
Thomas Witkowski committed
23
24
25
  {
    void *ctx;
    PCShellGetContext(pc, &ctx);
26
    NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
Thomas Witkowski's avatar
Thomas Witkowski committed
27

28
29
30
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
31
  }
32
  
33
34
35
36
37

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
      massMatrixSolver(NULL),
38
39
40
41
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
      solution(NULL)
42
43
44
45
46
47
48
49
50
51
52
  {
    Parameters::get(initFileStr + "->stokes->pressure component", 
		    pressureComponent);

    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");

    TEST_EXIT(pressureComponent == 2)("Fixed for pressure component = 2\n");
  }


53
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
54
55
56
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

57
    // Create FGMRES based outer solver
58
    KSPCreate(mpiCommGlobal, &ksp);
59
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
60
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
61
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
62

63
64
    // Create null space information.
    setConstantNullSpace(ksp, pressureComponent);
65
66
67
68
69
70
71
  }


  void PetscSolverNavierStokes::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");

72
73
74
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
75

76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);

    KSPSetUp(kspInterior);

    KSP *subKsp;
    int nSubKsp;
    PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);

    TEST_EXIT(nSubKsp == 2)
      ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");

    KSP kspVelocity = subKsp[0];
    KSP kspSchur = subKsp[1];
    PetscFree(subKsp);

    petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);

    KSPSetType(kspSchur, KSPPREONLY);
98
    PC pcSub;
99
100
101
102
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
103

104
105
106
107
    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
    KSPSetNullSpace(kspSchur, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
108

109
110

    // === Mass matrix solver ===
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
    Simple_ZOT zot;
    massOp.addTerm(&zot);
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
    Simple_SOT sot;
    laplaceOp.addTerm(&sot);
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


    // === Create convection-diffusion operator ===

135
136
137
138
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
139
140
141

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
142
    Simple_ZOT conDif0(*invTau);
143
    conDifOp.addTerm(&conDif0);
144
    Simple_SOT conDif1(*nu);
145
    conDifOp.addTerm(&conDif1);
146
    VecAtQP_FOT conDif2(&vx, &idFct, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
147
    conDifOp.addTerm(&conDif2, GRD_PHI);
148
    VecAtQP_FOT conDif3(&vy, &idFct, 1);
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    conDifOp.addTerm(&conDif3, GRD_PHI);
150
151
152
153
154
155

    conDifMatrix.assembleOperator(conDifOp);

    conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
    conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);

156
157

    // === Setup solver ===
158
159
160
161
162

    matShellContext.kspMass = massMatrixSolver->getSolver();
    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
    matShellContext.matConDif = conDifMatrixSolver->getMatInterior();    

163
164
    petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
    petsc_helper::setSolver(matShellContext.kspLaplace, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
165
166
  }

167

168
}