PetscSolverNavierStokes.cc 6.67 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
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
32
  }
33
  
34
35
36
37
38

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
      massMatrixSolver(NULL),
39
40
41
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
42
43
      solution(NULL),
      phase(NULL)
44
  {
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
45
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
46
47
48
49
50
51
52
53
54
		    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");
  }


55
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
56
57
58
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

59
    // Create FGMRES based outer solver
60
    KSPCreate(mpiCommGlobal, &ksp);
61
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
62
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
63
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
64
    
65
    // Create null space information.
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
66
    setConstantNullSpace(ksp, pressureComponent, true);
67
68
69
70
71
72
73
  }


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

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
74
75
76
77
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");

78
79
80
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
81

82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    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);

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
101
    petsc_helper::setSolver(kspVelocity, "velocity_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
102
103

    KSPSetType(kspSchur, KSPPREONLY);
104
    PC pcSub;
105
106
107
108
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
109

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
110
    setConstantNullSpace(kspSchur);
111

112
113

    // === Mass matrix solver ===
114
115
116
117

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
118
119
120
121
    //    if (!phase)
      massOp.addTerm(new Simple_ZOT);
    // else
    //   massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
122
123
124
125
126
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


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

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

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
133
134
135
136
    if (!phase)
      laplaceOp.addTerm(new Simple_SOT);
    else
      laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct));
137
138
139
140
141
142
143
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

144
145
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
146
    DOFVector<double> vp(pressureFeSpace, "vp");
147
148
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
149
150
151

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178

    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);
    }
179
180
181
182
183
184

    conDifMatrix.assembleOperator(conDifOp);

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

185
186

    // === Setup solver ===
187
188
189
190
191

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

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
    petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
    //    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);

    KSPSetType(matShellContext.kspLaplace, KSPRICHARDSON);
    KSPSetTolerances(matShellContext.kspLaplace, 0.0, 1e-14, PETSC_DEFAULT, 1);
    KSPSetOptionsPrefix(matShellContext.kspLaplace, "laplace_");
    KSPSetFromOptions(matShellContext.kspLaplace);

    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
    KSPSetNullSpace(matShellContext.kspLaplace, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);

    {
      PC pc;
      KSPGetPC(matShellContext.kspLaplace, &pc);
      PCSetType(pc, PCHYPRE);
      PCSetFromOptions(pc);
    }



    //    setConstantNullSpace(matShellContext.kspLaplace);

216
217
  }

218

219
}