PetscSolverNavierStokes.cc 11.7 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
28
#include "tools.h"

29
namespace AMDiS { namespace Parallel {
30
31
32

  using namespace std;

33

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

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

50
51
52
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
53
54
    
    PetscFunctionReturn(0);
55
  }
56
  
57
58

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

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

    Parameters::get(initFileStr + "->navierstokes->use old initial guess", 
		    useOldInitialGuess);
84
85
86
87
88
89
90
91
92

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

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

    Parameters::get(initFileStr + "->navierstokes->laplace solver", 
		    laplaceSolutionMode);
93
94
95
96
97
98
99
100
  }


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

101
    if (useOldInitialGuess) {      
102
103
104
105
106
107
108
109
110
111
      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);
112
113
114
  }


115
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
116
117
118
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

119
    printMem("begin(initSolver)"); 
120
    // Create FGMRES based outer solver
121
122
    
    MSG("CREATE POS 1: %p\n", &ksp);
123
    KSPCreate(domainComm, &ksp);
124
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
125
126
127
128
    if (getInfo() >= 10)
      KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
    else if (getInfo() >= 20)
      KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
129
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
130
    
131
    // Create null space information.
132
133
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
134
135
    
    printMem("end(initSolver)"); 
136
137
138
139
140
141
142
  }


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

143
    Timer t;
144

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
145
146
147
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");
148
149
150
    
    printMem("begin(initPreconditioner)"); 
    
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
151

Thomas Witkowski's avatar
Thomas Witkowski committed
152
153
    int dim = componentSpaces[pressureComponent]->getMesh()->getDim();

154
    vector<int> velocityComponents;
155
    for (int i = 0; i < dim; i++)
156
      velocityComponents.push_back(i);
157

158
159
160
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
161

162
163
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);
164
    PCSetFromOptions(pc);
165

Thomas Witkowski's avatar
Thomas Witkowski committed
166
    KSPSetUp(kspInterior);
167
168
169
170
171

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

172

173
174
175
176
177
178
179
    TEST_EXIT(nSubKsp == 2)
      ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");

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

180
181
182
    switch (velocitySolutionMode) {
    case 0:      
      petsc_helper::setSolver(kspVelocity, "", 
183
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
184
185
186
187
188
189
190
191
192
      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);
    }
    
193
194

    KSPSetType(kspSchur, KSPPREONLY);
195
    PC pcSub;
196
197
198
199
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
200

201
202
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
203

204
    // === Mass matrix solver ===
205
206
207

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
208
209
    {
      Operator massOp(pressureFeSpace, pressureFeSpace);
210
      ZeroOrderTerm *massTerm = nullptr;
211
      if ((!phase) || (*nu == 0.0))
212
213
	massTerm = new Simple_ZOT;
      else
214
        massTerm = new VecAtQP_ZOT(phase);
215
216
217
218
      massOp.addTerm(massTerm);
      massMatrix.assembleOperator(massOp);
      delete massTerm;
    }
219
220
221
222
223
224
225
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
226
227
    {
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
228
      SecondOrderTerm *laplaceTerm = nullptr;      
229
      if ((!phase) || (*nu == 0.0))
230
231
	laplaceTerm = new Simple_SOT;
      else
232
	laplaceTerm = new VecAtQP_SOT(phase);
233
234
235
236
      laplaceOp.addTerm(laplaceTerm);
      laplaceMatrix.assembleOperator(laplaceOp);
      delete laplaceTerm;
    }
237
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
238
239
240
241
242
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

243
244
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
Thomas Witkowski committed
245
    DOFVector<double> vz(pressureFeSpace, "vz");
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
246
    DOFVector<double> vp(pressureFeSpace, "vp");
247
    vx.interpol(solution->getDOFVector(0));
248
249
250
    if (dim >= 2)
      vy.interpol(solution->getDOFVector(1));
    if (dim >= 3)
Thomas Witkowski's avatar
Thomas Witkowski committed
251
      vz.interpol(solution->getDOFVector(2));
252
253

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

255
256
257
    {
      Operator conDifOp(pressureFeSpace, pressureFeSpace);
      
258
259
260
      ZeroOrderTerm *conDif0 = nullptr;
      SecondOrderTerm *conDif1 = nullptr;
      FirstOrderTerm *conDif2 = nullptr, *conDif3 = nullptr, *conDif4 = nullptr;
261
262
263
264
265
266
267
268

      if (!phase) {
	MSG("INIT WITHOUT PHASE!\n");
	
	conDif0 = new Simple_ZOT(*invTau);
	conDifOp.addTerm(conDif0);
	conDif1 = new Simple_SOT(*nu);
	conDifOp.addTerm(conDif1);
269
	conDif2 = new VecAtQP_FOT(&vx, nullptr, 0);
270
	conDifOp.addTerm(conDif2, GRD_PHI);
271
	if (dim >= 2) {
272
	  conDif3 = new VecAtQP_FOT(&vy, nullptr, 1);
273
274
	  conDifOp.addTerm(conDif3, GRD_PHI);
	}
275
	if (dim == 3) {
276
	  conDif4 = new VecAtQP_FOT(&vz, nullptr, 2);	
277
278
	  conDifOp.addTerm(conDif4, GRD_PHI);
	}
279
      } else { // no phase given
280
281
282
	
	vp.interpol(phase);
	
283
	if (*nu > 0.0) {
284
	  conDif0 = new VecAtQP_ZOT(&vp, nullptr, *invTau);
285
	  conDifOp.addTerm(conDif0);	
286
	  conDif1 = new VecAtQP_SOT(&vp, nullptr, *nu);
287
	  conDifOp.addTerm(conDif1);	
288
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, nullptr, 0);
289
290
	  conDifOp.addTerm(conDif2, GRD_PHI);
	  
291
	  if (dim >= 2) {
292
	    conDif3 = new Vec2AtQP_FOT(&vy, &vp, nullptr, 1);
293
294
	    conDifOp.addTerm(conDif3, GRD_PHI);
	  }
295
	  if (dim == 3) {
296
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, nullptr, 2);
297
298
299
300
301
302
303
304
305
306
	    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);
	  
307
308
309
310
	  if (dim >= 2) {
	    conDif3 = new Vec2AtQP_FOT(&vy, &vp, new LinearInterpolation2(*rho1,*rho2), 1);
	    conDifOp.addTerm(conDif3, GRD_PHI);
	  }
311
312
313
314
	  if (dim == 3) {
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, new LinearInterpolation2(*rho1,*rho2), 2);
	    conDifOp.addTerm(conDif4, GRD_PHI);
	  }
315
	}
316
      }/**/
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
317

318
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
319

320
321
322
      delete conDif0;
      delete conDif1;
      delete conDif2;
323
324
      if (dim >= 2)
	delete conDif3;
325
326
      if (dim == 3)
	delete conDif4;
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
327
    }
328
329
330
331

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

332
333

    // === Setup solver ===
334
335
336
337
338

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

339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
    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_", 
355
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
356
357
      break;
    case 1:
358
      petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "laplace_", 
359
360
361
362
363
364
				    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
365

366
    setConstantNullSpace(matShellContext.kspLaplace);
367
368
369
    
    
    printMem("end(initPreconditioner)"); 
370
371

    MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n", 
372
	t.elapsed());
373
374
  }

375

376
377
378
379
380
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
381
    massMatrixSolver->destroyVectorData();
382
    laplaceMatrixSolver->destroyMatrixData();
383
    laplaceMatrixSolver->destroyVectorData();
384
    conDifMatrixSolver->destroyMatrixData();
385
    conDifMatrixSolver->destroyVectorData();
386

387
388
389
390
391
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    conDifMatrixSolver->destroyVectorData();
    
    
392
    delete massMatrixSolver;
393
    massMatrixSolver = nullptr;
394
395

    delete laplaceMatrixSolver;
396
    laplaceMatrixSolver = nullptr;
397
398

    delete conDifMatrixSolver;
399
    conDifMatrixSolver = nullptr;
400
  }
401
} }