PetscSolverNavierStokes.cc 11.2 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * 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.
18
 *
19
 ******************************************************************************/
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
    // Project out constant null space
    {
      int vecSize;
      VecGetSize(x, &vecSize);
      PetscScalar vecSum;
      VecSum(x, &vecSum);
      vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
45
      VecShift(x, vecSum);
46 47
    }

48 49 50
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
51

52
    PetscFunctionReturn(0);
53
  }
54

55 56

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

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

80
    Parameters::get(initFileStr + "->navierstokes->use old initial guess",
81
		    useOldInitialGuess);
82

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

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

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


94
  void PetscSolverNavierStokes::solvePetscMatrix(SystemVector &vec,
95 96 97 98
						 AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverNavierStokes::solvePetscMatrix()");

99
    if (useOldInitialGuess) {
100
      VecSet(getVecSolInterior(), 0.0);
101

102 103
      for (int i = 0; i < solution->getSize(); i++)
	setDofVector(getVecSolInterior(), solution->getDOFVector(i), i, true);
104

105 106 107 108 109
      vecSolAssembly();
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
    }

    PetscSolverGlobalMatrix::solvePetscMatrix(vec, adaptInfo);
110 111 112
  }


113
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
114 115
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");
116
    // Create FGMRES based outer solver
117

118
    MSG("CREATE POS 1: %p\n", &ksp);
119
    KSPCreate(domainComm, &ksp);
120
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
121
    if (getInfo() >= 10)
122
      petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
123
    else if (getInfo() >= 20)
124
      petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
125
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
126

127
    // Create null space information.
128 129
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
130
    }
131 132 133 134 135 136


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

137
    Timer t;
138

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
139 140
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
141
    TEST_EXIT(solution)("solution pointer not set!\n");
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
142

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

145
    vector<int> velocityComponents;
146
    for (int i = 0; i < dim; i++)
147
      velocityComponents.push_back(i);
148

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

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

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

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

163

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

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

171
    switch (velocitySolutionMode) {
172 173
    case 0:
      petsc_helper::setSolver(kspVelocity, "",
174
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
175 176
      break;
    case 1:
177
      petsc_helper::setSolverWithLu(kspVelocity, "", KSPPREONLY,
178 179 180 181 182
				    PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
    }
183

184 185

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

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

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

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


    // === Laplace matrix solver ===

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


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

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

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

246 247
    {
      Operator conDifOp(pressureFeSpace, pressureFeSpace);
248

249 250 251
      ZeroOrderTerm *conDif0 = NULL;
      SecondOrderTerm *conDif1 = NULL;
      FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;
252 253 254

      if (!phase) {
	MSG("INIT WITHOUT PHASE!\n");
255

256 257 258 259
	conDif0 = new Simple_ZOT(*invTau);
	conDifOp.addTerm(conDif0);
	conDif1 = new Simple_SOT(*nu);
	conDifOp.addTerm(conDif1);
260
	conDif2 = new VecAtQP_FOT(&vx, NULL, 0);
261
	conDifOp.addTerm(conDif2, GRD_PHI);
262
	if (dim >= 2) {
263
	  conDif3 = new VecAtQP_FOT(&vy, NULL, 1);
264 265
	  conDifOp.addTerm(conDif3, GRD_PHI);
	}
266
	if (dim == 3) {
267
	  conDif4 = new VecAtQP_FOT(&vz, NULL, 2);
268 269
	  conDifOp.addTerm(conDif4, GRD_PHI);
	}
270
      } else { // no phase given
271

272
	vp.interpol(phase);
273

274
	if (*nu > 0.0) {
275
	  conDif0 = new VecAtQP_ZOT(&vp, NULL, *invTau);
276
	  conDifOp.addTerm(conDif0);
277
	  conDif1 = new VecAtQP_SOT(&vp, NULL, *nu);
278
	  conDifOp.addTerm(conDif1);
279
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, NULL, 0);
280
	  conDifOp.addTerm(conDif2, GRD_PHI);
281

282
	  if (dim >= 2) {
283
	    conDif3 = new Vec2AtQP_FOT(&vy, &vp, NULL, 1);
284 285
	    conDifOp.addTerm(conDif3, GRD_PHI);
	  }
286
	  if (dim == 3) {
287
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, NULL, 2);
288 289
	    conDifOp.addTerm(conDif4, GRD_PHI);
	  }
290
	} else {
291
	  conDif0 = new VecAtQP_ZOT(&vp, new LinearInterpolation(*rho1,*rho2,*invTau));
292
	  conDifOp.addTerm(conDif0);
293
	  conDif1 = new VecAtQP_SOT(&vp, new LinearInterpolation(*nu1,*nu2));
294
	  conDifOp.addTerm(conDif1);
295 296
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, new LinearInterpolation2(*rho1,*rho2), 0);
	  conDifOp.addTerm(conDif2, GRD_PHI);
297

298 299 300 301
	  if (dim >= 2) {
	    conDif3 = new Vec2AtQP_FOT(&vy, &vp, new LinearInterpolation2(*rho1,*rho2), 1);
	    conDifOp.addTerm(conDif3, GRD_PHI);
	  }
302 303 304 305
	  if (dim == 3) {
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, new LinearInterpolation2(*rho1,*rho2), 2);
	    conDifOp.addTerm(conDif4, GRD_PHI);
	  }
306
	}
307
      }/**/
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
308

309
      conDifMatrix.assembleOperator(conDifOp);
Thomas Witkowski's avatar
Thomas Witkowski committed
310

311 312 313
      delete conDif0;
      delete conDif1;
      delete conDif2;
314 315
      if (dim >= 2)
	delete conDif3;
316 317
      if (dim == 3)
	delete conDif4;
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
318
    }
319 320 321 322

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

323 324

    // === Setup solver ===
325 326 327

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

330 331
    switch (massSolutionMode) {
    case 0:
332
      petsc_helper::setSolver(matShellContext.kspMass, "mass_",
333 334 335
			      KSPCG, PCJACOBI, 0.0, 1e-14, 2);
      break;
    case 1:
336
      petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON,
337 338 339 340 341 342 343 344
				    PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode);
    }

    switch (laplaceSolutionMode) {
    case 0:
345
      petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_",
346
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
347 348
      break;
    case 1:
349 350
      petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "laplace_",
				    KSPRICHARDSON, PCLU, MATSOLVERMUMPS,
351 352 353 354 355
				    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
356

357
    setConstantNullSpace(matShellContext.kspLaplace);
358 359

    MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n",
360
	t.elapsed());
361 362
  }

363

364 365 366 367 368
  void PetscSolverNavierStokes::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
369
    massMatrixSolver->destroyVectorData();
370
    laplaceMatrixSolver->destroyMatrixData();
371
    laplaceMatrixSolver->destroyVectorData();
372
    conDifMatrixSolver->destroyMatrixData();
373
    conDifMatrixSolver->destroyVectorData();
374

375 376 377
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    conDifMatrixSolver->destroyVectorData();
378 379


380
    delete massMatrixSolver;
381
    massMatrixSolver = NULL;
382 383

    delete laplaceMatrixSolver;
384
    laplaceMatrixSolver = NULL;
385 386

    delete conDifMatrixSolver;
387
    conDifMatrixSolver = NULL;
388
  }
389
} }