PetscSolverNavierStokes.cc 9.84 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
108
    
    MSG("CREATE POS 1: %p\n", &ksp);
109
    KSPCreate(domainComm, &ksp);
110
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
111
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
112
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
113
    
114
    // Create null space information.
115
116
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
117
118
119
120
121
122
123
  }


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

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

130
131
132
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
    if (dim == 3)
      velocityComponents.push_back(2);
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);
141
    PCSetFromOptions(pc);
142
143
144
145
146
147
148
149
150
151
152
153
154
155

    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);

156
157
158
159
160
161
162
163
164
165
166
167
168
    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);
    }
    
169
170

    KSPSetType(kspSchur, KSPPREONLY);
171
    PC pcSub;
172
173
174
175
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
176

177
178
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
179

180
    // === Mass matrix solver ===
181
182
183

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
184
185
186
187
188
189
190
191
192
193
194
    {
      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;
    }
195
196
197
198
199
200
201
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
202
203
204
205
206
207
208
209
210
211
212
    {
      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;
    }
213
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
214
215
216
217
218
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

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

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
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
270
271
    {
      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
272

273
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
274

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

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

286
287

    // === Setup solver ===
288
289
290
291
292

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

293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    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
319

320
    setConstantNullSpace(matShellContext.kspLaplace);
321
322
  }

323

324
325
326
327
328
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
329
    massMatrixSolver->destroyVectorData();
330
    laplaceMatrixSolver->destroyMatrixData();
331
    laplaceMatrixSolver->destroyVectorData();
332
    conDifMatrixSolver->destroyMatrixData();
333
    conDifMatrixSolver->destroyVectorData();
334
335
336
337
338
339
340
341
342
343

    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
344
}