PetscSolverNavierStokes.cc 10.9 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
  }


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

88
    if (useOldInitialGuess) {      
89
90
91
92
93
94
95
96
97
98
      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

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    KSPSetUp(kspInterior);
148
149
150
151
152

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

153

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

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

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

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

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

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

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
189
190
191
    {
      Operator massOp(pressureFeSpace, pressureFeSpace);
      ZeroOrderTerm *massTerm = NULL;
192
      if ((!phase)||(*nu==0.0))
193
194
	massTerm = new Simple_ZOT;
      else
195
        massTerm = new VecAtQP_ZOT(phase, &idFct);
196
197
198
199
      massOp.addTerm(massTerm);
      massMatrix.assembleOperator(massOp);
      delete massTerm;
    }
200
201
202
203
204
205
206
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
207
208
209
    {
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
      SecondOrderTerm *laplaceTerm = NULL;      
210
      if ((!phase)||(*nu==0.0))
211
212
213
214
215
216
217
	laplaceTerm = new Simple_SOT;
      else
	laplaceTerm = new VecAtQP_SOT(phase, &idFct);
      laplaceOp.addTerm(laplaceTerm);
      laplaceMatrix.assembleOperator(laplaceOp);
      delete laplaceTerm;
    }
218
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
219
220
221
222
223
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

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

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

235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
    {
      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);
	}
258
259
        } else {
	MSG("INIT WITH PHASE: %e %e %e %e %e %e \n", *nu, *invTau, *rho1, *rho2, *nu1, *nu2);
260
261
262
	
	vp.interpol(phase);
	
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
	if (*nu>0.0) {
	  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);
	  }
	} else {	  
	  conDif0 = new VecAtQP_ZOT(&vp, new LinearInterpolation(*rho1,*rho2,*invTau));
	  conDifOp.addTerm(conDif0);	
	  conDif1 = new VecAtQP_SOT(&vp, new LinearInterpolation(*nu1,*nu2));
	  conDifOp.addTerm(conDif1);	
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, new LinearInterpolation2(*rho1,*rho2), 0);
	  conDifOp.addTerm(conDif2, GRD_PHI);
	  
	  conDif3 = new Vec2AtQP_FOT(&vy, &vp, new LinearInterpolation2(*rho1,*rho2), 1);
	  conDifOp.addTerm(conDif3, GRD_PHI);
	  
	  if (dim == 3) {
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, new LinearInterpolation2(*rho1,*rho2), 2);
	    conDifOp.addTerm(conDif4, GRD_PHI);
	  }
293
	}
294
      }/**/
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
295

296
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
297

298
299
300
301
302
303
      delete conDif0;
      delete conDif1;
      delete conDif2;
      delete conDif3;
      if (dim == 3)
	delete conDif4;
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
304
    }
305
306
307
308

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

309
310

    // === Setup solver ===
311
312
313
314
315

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

316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
    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_", 
332
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
333
334
335
336
337
338
339
340
341
      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
342

343
    setConstantNullSpace(matShellContext.kspLaplace);
344
345
346

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

349

350
351
352
353
354
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
355
    massMatrixSolver->destroyVectorData();
356
    laplaceMatrixSolver->destroyMatrixData();
357
    laplaceMatrixSolver->destroyVectorData();
358
    conDifMatrixSolver->destroyMatrixData();
359
    conDifMatrixSolver->destroyVectorData();
360

361
362
363
364
365
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    conDifMatrixSolver->destroyVectorData();
    
    
366
367
368
369
370
371
372
373
374
    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
375
}