PetscSolverNavierStokes.cc 9.65 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
51
52
      velocitySolutionMode(0),
      massSolutionMode(0),
      laplaceSolutionMode(0),
53
      massMatrixSolver(NULL),
54
55
56
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
57
58
      solution(NULL),
      phase(NULL)
59
  {
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
60
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
61
62
63
		    pressureComponent);
    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");
64
65
66
67

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

    Parameters::get(initFileStr + "->navierstokes->use old initial guess", 
		    useOldInitialGuess);
71
72
73
74
75
76
77
78
79

    Parameters::get(initFileStr + "->navierstokes->velocity solver", 
		    velocitySolutionMode);

    Parameters::get(initFileStr + "->navierstokes->mass solver", 
		    massSolutionMode);

    Parameters::get(initFileStr + "->navierstokes->laplace solver", 
		    laplaceSolutionMode);
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
  }


  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);
99
100
101
  }


102
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
103
104
105
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

106
    // Create FGMRES based outer solver
107
    KSPCreate(mpiCommGlobal, &ksp);
108
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
109
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
110
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
111
    
112
    // Create null space information.
113
114
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
115
116
117
118
119
120
121
  }


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

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
122
123
124
125
    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
126
127
    int dim = componentSpaces[pressureComponent]->getMesh()->getDim();

128
129
130
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
    if (dim == 3)
      velocityComponents.push_back(2);
133

134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    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);

153
154
155
156
157
158
159
160
161
162
163
164
165
    switch (velocitySolutionMode) {
    case 0:      
      petsc_helper::setSolver(kspVelocity, "", 
			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
      break;
    case 1:
      petsc_helper::setSolverWithLu(kspVelocity, "", KSPPREONLY, 
				    PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
    }
    
166
167

    KSPSetType(kspSchur, KSPPREONLY);
168
    PC pcSub;
169
170
171
172
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
173

174
175
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
176

177
178

    // === Mass matrix solver ===
179
180
181

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
182
183
184
185
186
187
188
189
190
191
192
    {
      Operator massOp(pressureFeSpace, pressureFeSpace);
      ZeroOrderTerm *massTerm = NULL;
      if (!phase)
	massTerm = new Simple_ZOT;
      else
	massTerm = new VecAtQP_ZOT(phase, &idFct);
      massOp.addTerm(massTerm);
      massMatrix.assembleOperator(massOp);
      delete massTerm;
    }
193
194
195
196
197
198
199
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
200
201
202
203
204
205
206
207
208
209
210
    {
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
      SecondOrderTerm *laplaceTerm = NULL;      
      if (!phase)
	laplaceTerm = new Simple_SOT;
      else
	laplaceTerm = new VecAtQP_SOT(phase, &idFct);
      laplaceOp.addTerm(laplaceTerm);
      laplaceMatrix.assembleOperator(laplaceOp);
      delete laplaceTerm;
    }
211
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
212
213
214
215
216
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

217
218
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    DOFVector<double> vz(pressureFeSpace, "vz");
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
220
    DOFVector<double> vp(pressureFeSpace, "vp");
221
222
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
    if (dim == 3)
      vz.interpol(solution->getDOFVector(2));
225
226

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
227

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
    {
      Operator conDifOp(pressureFeSpace, pressureFeSpace);
      
      ZeroOrderTerm *conDif0 = NULL;
      SecondOrderTerm *conDif1 = NULL;
      FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;

      if (!phase) {
	MSG("INIT WITHOUT PHASE!\n");
	
	conDif0 = new Simple_ZOT(*invTau);
	conDifOp.addTerm(conDif0);
	conDif1 = new Simple_SOT(*nu);
	conDifOp.addTerm(conDif1);
	conDif2 = new VecAtQP_FOT(&vx, &idFct, 0);
	conDifOp.addTerm(conDif2, GRD_PHI);
	conDif3 = new VecAtQP_FOT(&vy, &idFct, 1);
	conDifOp.addTerm(conDif3, GRD_PHI);
	
	if (dim == 3) {
	  conDif4 = new VecAtQP_FOT(&vz, &idFct, 2);	
	  conDifOp.addTerm(conDif4, GRD_PHI);
	}
      } else {
	MSG("INIT WITH PHASE: %e %e\n", *nu, *invTau);
	
	vp.interpol(phase);
	conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau));
	conDifOp.addTerm(conDif0);	
	conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu));
	conDifOp.addTerm(conDif1);	
	conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0);
	conDifOp.addTerm(conDif2, GRD_PHI);
	
	conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1);
	conDifOp.addTerm(conDif3, GRD_PHI);
	
	if (dim == 3) {
	  conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2);
	  conDifOp.addTerm(conDif4, GRD_PHI);
	}
      }
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
270

271
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
272

273
274
275
276
277
278
      delete conDif0;
      delete conDif1;
      delete conDif2;
      delete conDif3;
      if (dim == 3)
	delete conDif4;
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
279
    }
280
281
282
283

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

284
285

    // === Setup solver ===
286
287
288
289
290

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

291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
    switch (massSolutionMode) {
    case 0:
      petsc_helper::setSolver(matShellContext.kspMass, "mass_", 
			      KSPCG, PCJACOBI, 0.0, 1e-14, 2);
      break;
    case 1:
      petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON, 
				    PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode);
    }

    switch (laplaceSolutionMode) {
    case 0:
      petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", 
			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
      break;
    case 1:
      petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "mass_", 
				    KSPRICHARDSON, PCLU, MATSOLVERMUMPS, 
				    0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode);
    }
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
317

318
    setConstantNullSpace(matShellContext.kspLaplace);
319
320
  }

321

322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
    laplaceMatrixSolver->destroyMatrixData();
    conDifMatrixSolver->destroyMatrixData();

    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
339
}