PetscSolverNavierStokes.cc 6.55 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"
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
16
#include "TransformDOF.h"
17
18
19
20
21

namespace AMDiS {

  using namespace std;

22

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

29
30
31
32
33
34
35
36
37
38
    // Project out constant null space
    {
      int vecSize;
      VecGetSize(x, &vecSize);
      PetscScalar vecSum;
      VecSum(x, &vecSum);
      vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
      VecShift(x, vecSum); 
    }

39
40
41
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
42
  }
43
  
44
45
46
47

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
48
      pressureNullSpace(true),
49
      massMatrixSolver(NULL),
50
51
52
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
53
54
      solution(NULL),
      phase(NULL)
55
  {
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
56
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
57
58
59
60
		    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");
61
62
63
64

    Parameters::get(initFileStr + "->navierstokes->pressure null space",
		    pressureNullSpace);
    TEST_EXIT(pressureNullSpace)("This is not yet tested, may be wrong!\n");
65
66
67
  }


68
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
69
70
71
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

72
    // Create FGMRES based outer solver
73
    KSPCreate(mpiCommGlobal, &ksp);
74
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
75
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
76
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
77
    
78
    // Create null space information.
79
80
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
81
82
83
84
85
86
87
  }


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

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
88
89
90
91
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");

92
93
94
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
95

96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
    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);

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

    KSPSetType(kspSchur, KSPPREONLY);
118
    PC pcSub;
119
120
121
122
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
123

124
125
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
126

127
128

    // === Mass matrix solver ===
129
130
131
132

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
133
    if (!phase)
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
134
      massOp.addTerm(new Simple_ZOT);
135
136
    else
      massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
137
138
139
140
141
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
142
143
    //    VtkWriter::writeFile(phase, "phase.vtu");

144
145
146
147
    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
148
149
150
151
    if (!phase)
      laplaceOp.addTerm(new Simple_SOT);
    else
      laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct));
152
153
154
155
156
157
158
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

159
160
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
161
    DOFVector<double> vp(pressureFeSpace, "vp");
162
163
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
164
165
166

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

    if (!phase) {
      MSG("INIT WITHOUT PHASE!\n");
      Simple_ZOT *conDif0 = new Simple_ZOT(*invTau);
      conDifOp.addTerm(conDif0);
      Simple_SOT *conDif1 = new Simple_SOT(*nu);
      conDifOp.addTerm(conDif1);
      VecAtQP_FOT *conDif2 = new VecAtQP_FOT(&vx, &idFct, 0);
      conDifOp.addTerm(conDif2, GRD_PHI);
      VecAtQP_FOT *conDif3 = new VecAtQP_FOT(&vy, &idFct, 1);
      conDifOp.addTerm(conDif3, GRD_PHI);
    } else {
      MSG("INIT WITH PHASE!\n");
      
      vp.interpol(phase);
      VecAtQP_ZOT *conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau));
      conDifOp.addTerm(conDif0);

      VecAtQP_SOT *conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu));
      conDifOp.addTerm(conDif1);

      Vec2AtQP_FOT *conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0);
      conDifOp.addTerm(conDif2, GRD_PHI);

      Vec2AtQP_FOT *conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1);
      conDifOp.addTerm(conDif3, GRD_PHI);
    }
194
195
196
197
198
199

    conDifMatrix.assembleOperator(conDifOp);

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

200
201

    // === Setup solver ===
202
203
204
205
206

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

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
207
    petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
208
    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
209

210
    setConstantNullSpace(matShellContext.kspLaplace);
211
212
  }

213

214
}