PetscSolverNavierStokes.cc 7.58 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
      useOldInitialGuess(false),
50
      massMatrixSolver(NULL),
51
52
53
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
54
55
      solution(NULL),
      phase(NULL)
56
  {
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
57
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
58
59
60
		    pressureComponent);
    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

    Parameters::get(initFileStr + "->navierstokes->use old initial guess", 
		    useOldInitialGuess);
  }


  void PetscSolverNavierStokes::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverNavierStokes::solvePetscMatrix()");

    if (useOldInitialGuess) {
      VecSet(getVecSolInterior(), 0.0);
      
      for (int i = 0; i < solution->getSize(); i++)
	setDofVector(getVecSolInterior(), solution->getDOFVector(i), i, true);
      
      vecSolAssembly();
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
    }

    PetscSolverGlobalMatrix::solvePetscMatrix(vec, adaptInfo);
87
88
89
  }


90
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
91
92
93
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

94
    // Create FGMRES based outer solver
95
    KSPCreate(mpiCommGlobal, &ksp);
96
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
97
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
98
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
99
    
100
    // Create null space information.
101
102
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
103
104
105
106
107
108
109
  }


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

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
110
111
112
113
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
    int dim = componentSpaces[pressureComponent]->getMesh()->getDim();

116
117
118
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
Thomas Witkowski's avatar
Thomas Witkowski committed
119
120
    if (dim == 3)
      velocityComponents.push_back(2);
121

122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
    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);

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

    KSPSetType(kspSchur, KSPPREONLY);
144
    PC pcSub;
145
146
147
148
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
149

150
151
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
152

153
154

    // === Mass matrix solver ===
155
156
157
158

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
159
    if (!phase)
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
160
      massOp.addTerm(new Simple_ZOT);
161
162
    else
      massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
163
164
165
166
167
168
169
170
171
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
172
173
174
175
    if (!phase)
      laplaceOp.addTerm(new Simple_SOT);
    else
      laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct));
176
    laplaceMatrix.assembleOperator(laplaceOp);
177
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
178
179
180
181
182
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

183
184
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
Thomas Witkowski committed
185
    DOFVector<double> vz(pressureFeSpace, "vz");
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
186
    DOFVector<double> vp(pressureFeSpace, "vp");
187
188
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
    if (dim == 3)
      vz.interpol(solution->getDOFVector(2));
191
192
193

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
194
195
196
197
198
199
200
201
202
203
204

    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);
Thomas Witkowski's avatar
Thomas Witkowski committed
205
206
207
208
209

      if (dim == 3) {
	VecAtQP_FOT *conDif4 = new VecAtQP_FOT(&vz, &idFct, 2);	
	conDifOp.addTerm(conDif4, GRD_PHI);
      }
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
    } 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);
Thomas Witkowski's avatar
Thomas Witkowski committed
225
226
227
228
229

      if (dim == 3) {
	Vec2AtQP_FOT *conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2);
	conDifOp.addTerm(conDif4, GRD_PHI);
      }
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
230
    }
231
232
233
234
235
236

    conDifMatrix.assembleOperator(conDifOp);

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

237
238

    // === Setup solver ===
239
240
241
242
243

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

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

247
    setConstantNullSpace(matShellContext.kspLaplace);
248
249
  }

250

251
}