PetscSolverNavierStokes.cc 9.97 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()");

124
125
126
    MPI::COMM_WORLD.Barrier();
    double wtime = MPI::Wtime();

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

133
134
135
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
    if (dim == 3)
      velocityComponents.push_back(2);
138

139
140
141
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
142

143
144
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);
145
    PCSetFromOptions(pc);
146
147
148
149
150
151


    KSP *subKsp;
    int nSubKsp;
    PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);

152

153
154
155
156
157
158
159
    TEST_EXIT(nSubKsp == 2)
      ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");

    KSP kspVelocity = subKsp[0];
    KSP kspSchur = subKsp[1];
    PetscFree(subKsp);

160
161
162
    switch (velocitySolutionMode) {
    case 0:      
      petsc_helper::setSolver(kspVelocity, "", 
163
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
164
165
166
167
168
169
170
171
172
      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);
    }
    
173
174

    KSPSetType(kspSchur, KSPPREONLY);
175
    PC pcSub;
176
177
178
179
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
180

181
182
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
183

184
    // === Mass matrix solver ===
185
186
187

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


    // === Laplace matrix solver ===

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


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

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

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

277
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
278

279
280
281
282
283
284
      delete conDif0;
      delete conDif1;
      delete conDif2;
      delete conDif3;
      if (dim == 3)
	delete conDif4;
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
285
    }
286
287
288
289

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

290
291

    // === Setup solver ===
292
293
294
295
296

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

297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
    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_", 
313
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
314
315
316
317
318
319
320
321
322
      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
323

324
    setConstantNullSpace(matShellContext.kspLaplace);
325
326
327

    MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n", 
	MPI::Wtime() - wtime);
328
329
  }

330

331
332
333
334
335
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
336
    massMatrixSolver->destroyVectorData();
337
    laplaceMatrixSolver->destroyMatrixData();
338
    laplaceMatrixSolver->destroyVectorData();
339
    conDifMatrixSolver->destroyMatrixData();
340
    conDifMatrixSolver->destroyVectorData();
341
342
343
344
345
346
347
348
349
350

    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
351
}