PetscSolverNavierStokes.cc 11.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21
22
23



#include "parallel/PetscSolverNavierStokes.h"
24
#include "parallel/PetscHelper.h"
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
25
#include "TransformDOF.h"
26

27
namespace AMDiS { namespace Parallel {
28
29
30

  using namespace std;

31

32
  PetscErrorCode pcSchurShell(PC pc, Vec x, Vec y)
Thomas Witkowski's avatar
Thomas Witkowski committed
33
34
35
  {
    void *ctx;
    PCShellGetContext(pc, &ctx);
36
    NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
Thomas Witkowski's avatar
Thomas Witkowski committed
37

38
39
40
41
42
43
44
45
46
47
    // 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); 
    }

48
49
50
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
51
  }
52
  
53
54
55
56

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
57
      pressureNullSpace(true),
58
      useOldInitialGuess(false),
59
60
61
      velocitySolutionMode(0),
      massSolutionMode(0),
      laplaceSolutionMode(0),
62
      massMatrixSolver(NULL),
63
64
65
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
66
67
      solution(NULL),
      phase(NULL)
68
  {
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
69
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
70
71
72
		    pressureComponent);
    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");
73
74
75
76

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

    Parameters::get(initFileStr + "->navierstokes->use old initial guess", 
		    useOldInitialGuess);
80
81
82
83
84
85
86
87
88

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

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

    Parameters::get(initFileStr + "->navierstokes->laplace solver", 
		    laplaceSolutionMode);
89
90
91
92
93
94
95
96
  }


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

97
    if (useOldInitialGuess) {      
98
99
100
101
102
103
104
105
106
107
      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);
108
109
110
  }


111
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
112
113
114
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

115
    // Create FGMRES based outer solver
116
117
    
    MSG("CREATE POS 1: %p\n", &ksp);
118
    KSPCreate(domainComm, &ksp);
119
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
120
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
121
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
122
    
123
    // Create null space information.
124
125
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
126
127
128
129
130
131
132
  }


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

133
134
135
    MPI::COMM_WORLD.Barrier();
    double wtime = MPI::Wtime();

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
136
137
138
139
    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
140
141
    int dim = componentSpaces[pressureComponent]->getMesh()->getDim();

142
143
144
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
Thomas Witkowski's avatar
Thomas Witkowski committed
145
146
    if (dim == 3)
      velocityComponents.push_back(2);
147

148
149
150
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
151

152
153
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);
154
    PCSetFromOptions(pc);
155

Thomas Witkowski's avatar
Thomas Witkowski committed
156
    KSPSetUp(kspInterior);
157
158
159
160
161

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

162

163
164
165
166
167
168
169
    TEST_EXIT(nSubKsp == 2)
      ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");

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

170
171
172
    switch (velocitySolutionMode) {
    case 0:      
      petsc_helper::setSolver(kspVelocity, "", 
173
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
174
175
176
177
178
179
180
181
182
      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);
    }
    
183
184

    KSPSetType(kspSchur, KSPPREONLY);
185
    PC pcSub;
186
187
188
189
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
190

191
192
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
193

194
    // === Mass matrix solver ===
195
196
197

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
198
199
200
    {
      Operator massOp(pressureFeSpace, pressureFeSpace);
      ZeroOrderTerm *massTerm = NULL;
201
      if ((!phase)||(*nu==0.0))
202
203
	massTerm = new Simple_ZOT;
      else
204
        massTerm = new VecAtQP_ZOT(phase, &idFct);
205
206
207
208
      massOp.addTerm(massTerm);
      massMatrix.assembleOperator(massOp);
      delete massTerm;
    }
209
210
211
212
213
214
215
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
216
217
218
    {
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
      SecondOrderTerm *laplaceTerm = NULL;      
219
      if ((!phase)||(*nu==0.0))
220
221
222
223
224
225
226
	laplaceTerm = new Simple_SOT;
      else
	laplaceTerm = new VecAtQP_SOT(phase, &idFct);
      laplaceOp.addTerm(laplaceTerm);
      laplaceMatrix.assembleOperator(laplaceOp);
      delete laplaceTerm;
    }
227
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
228
229
230
231
232
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

233
234
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    DOFVector<double> vz(pressureFeSpace, "vz");
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
236
    DOFVector<double> vp(pressureFeSpace, "vp");
237
238
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
    if (dim == 3)
      vz.interpol(solution->getDOFVector(2));
241
242

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

244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
    {
      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);
	}
267
        } else {
268
269
270
	
	vp.interpol(phase);
	
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
	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);
	  }
301
	}
302
      }/**/
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
303

304
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
305

306
307
308
309
310
311
      delete conDif0;
      delete conDif1;
      delete conDif2;
      delete conDif3;
      if (dim == 3)
	delete conDif4;
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
312
    }
313
314
315
316

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

317
318

    // === Setup solver ===
319
320
321
322
323

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

324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
    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_", 
340
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
341
342
343
344
345
346
347
348
349
      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
350

351
    setConstantNullSpace(matShellContext.kspLaplace);
352
353
354

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

357

358
359
360
361
362
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
363
    massMatrixSolver->destroyVectorData();
364
    laplaceMatrixSolver->destroyMatrixData();
365
    laplaceMatrixSolver->destroyVectorData();
366
    conDifMatrixSolver->destroyMatrixData();
367
    conDifMatrixSolver->destroyVectorData();
368

369
370
371
372
373
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    conDifMatrixSolver->destroyVectorData();
    
    
374
375
376
377
378
379
380
381
382
    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
383
} }