PetscSolverNavierStokes2.cc 12.4 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
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:
Praetorius, Simon's avatar
Praetorius, Simon committed
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
 *
Praetorius, Simon's avatar
Praetorius, Simon committed
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
 ******************************************************************************/



#include "preconditioner/PetscSolverNavierStokes2.h"
#include "parallel/PetscHelper.h"
#include "TransformDOF.h"
#include "DirichletBC.h"
#include "Functors.h"

namespace AMDiS { namespace Parallel {

  using namespace std;


  PetscErrorCode pcSchurShell2(PC pc, Vec x, Vec y)
  {
    void *ctx;
    PCShellGetContext(pc, &ctx);
    NavierStokesSchurData2* data = static_cast<NavierStokesSchurData2*>(ctx);

    // Project out constant null space
    {
      int vecSize;
      VecGetSize(x, &vecSize);
      PetscScalar vecSum;
      VecSum(x, &vecSum);
      vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
47
      VecShift(x, vecSum);
Praetorius, Simon's avatar
Praetorius, Simon committed
48 49 50 51 52
    }

    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
53

Praetorius, Simon's avatar
Praetorius, Simon committed
54 55
    PetscFunctionReturn(0);
  }
56

Praetorius, Simon's avatar
Praetorius, Simon committed
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72

  PetscSolverNavierStokes2::PetscSolverNavierStokes2(string name)
    : PetscSolverGlobalMatrix(name, false),
      pressureComponent(-1),
      pressureNullSpace(true),
      useOldInitialGuess(false),
      velocitySolutionMode(0),
      massSolutionMode(0),
      laplaceSolutionMode(0),
      massMatrixSolver(NULL),
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
      solution(NULL),
      phase(NULL)
  {
73
    Parameters::get(name + "->navierstokes->pressure component",
Praetorius, Simon's avatar
Praetorius, Simon committed
74 75 76 77 78 79 80 81
		    pressureComponent);
    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");

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

82
    Parameters::get(name + "->navierstokes->use old initial guess",
Praetorius, Simon's avatar
Praetorius, Simon committed
83 84
		    useOldInitialGuess);

85
    Parameters::get(name + "->navierstokes->velocity solver",
Praetorius, Simon's avatar
Praetorius, Simon committed
86 87
		    velocitySolutionMode);

88
    Parameters::get(name + "->navierstokes->mass solver",
Praetorius, Simon's avatar
Praetorius, Simon committed
89 90
		    massSolutionMode);

91
    Parameters::get(name + "->navierstokes->laplace solver",
Praetorius, Simon's avatar
Praetorius, Simon committed
92
		    laplaceSolutionMode);
93 94 95


    Parameters::get(name + "->navierstokes->neumann boundary indices",
Praetorius, Simon's avatar
Praetorius, Simon committed
96 97
		    neumannBC);
  }
98

Praetorius, Simon's avatar
Praetorius, Simon committed
99 100 101 102 103 104 105
   int PetscSolverNavierStokes2::solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
				    SystemVector& x,
				    SystemVector& b,
				    bool createMatrixData,
				    bool storeMatrixData)
  {
    FUNCNAME("PetscSolverNavierStokes2::solveLinearSystem()");
106

Praetorius, Simon's avatar
Praetorius, Simon committed
107 108 109 110 111 112 113 114 115 116 117
    TEST_EXIT(meshDistributor)("No meshDistributor provided. Should not happen!\n");

    MPI::COMM_WORLD.Barrier();
    Timer t;

    systemMat = A.getOriginalMat();
    if (createMatrixData)
      fillPetscMatrix(const_cast< Matrix< DOFMatrix* >* >(systemMat));

    fillPetscRhs(&b);

118
    INFO(info, 8)("creation of parallel data structures needed %.5f seconds\n",
Praetorius, Simon's avatar
Praetorius, Simon committed
119 120 121 122 123 124 125 126
		  t.elapsed());

    solvePetscMatrix(x, NULL);

    if (!storeMatrixData) {
      destroyVectorData();
      destroyMatrixData();
    }
127

Praetorius, Simon's avatar
Praetorius, Simon committed
128 129 130 131 132 133 134 135
    return 0;
  }


  void PetscSolverNavierStokes2::initSolver(KSP &ksp)
  {
    // Create FGMRES based outer solver
    KSPCreate(domainComm, &ksp);
136 137 138
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp, getMatInterior(), getMatInterior());
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
139
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
140
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
141
    if (getInfo() >= 10)
142
      KSPMonitorSet(ksp, PETSC_MONITOR_CAST(KSPMonitorDefault), PETSC_NULL, PETSC_NULL);
Praetorius, Simon's avatar
Praetorius, Simon committed
143
    else if (getInfo() >= 20)
144
      KSPMonitorSet(ksp, PETSC_MONITOR_CAST(KSPMonitorTrueResidualNorm), PETSC_NULL, PETSC_NULL);
Praetorius, Simon's avatar
Praetorius, Simon committed
145
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
146

Praetorius, Simon's avatar
Praetorius, Simon committed
147 148 149
    // Create null space information.
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
150

Praetorius, Simon's avatar
Praetorius, Simon committed
151 152 153
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }
154

Praetorius, Simon's avatar
Praetorius, Simon committed
155 156
  void PetscSolverNavierStokes2::addDirichletBC(DOFMatrix& m, BoundaryType b)
  {
157
    DirichletBC *dirichletApply =
Praetorius, Simon's avatar
Praetorius, Simon committed
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
      new DirichletBC(b, new AMDiS::Const<double, WorldVector<double> >(0.0), m.getRowFeSpace(), m.getColFeSpace(), true);

    m.getBoundaryManager()->addBoundaryCondition(dirichletApply);
  }


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

    Timer t;

    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");

    int dow = Global::getGeo(WORLD);

    vector<int> velocityComponents;
    for (size_t i = 0; i < static_cast<size_t>(dow); i++)
      velocityComponents.push_back(i);

    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);

    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);
    PCSetFromOptions(pc);

    KSPSetUp(kspInterior);

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


    TEST_EXIT(nSubKsp == 2)
      ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");

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

    switch (velocitySolutionMode) {
203 204
    case 0:
      petsc_helper::setSolver(kspVelocity, "",
Praetorius, Simon's avatar
Praetorius, Simon committed
205 206 207
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
      break;
    case 1:
208
      petsc_helper::setSolverWithLu(kspVelocity, "", KSPPREONLY,
Praetorius, Simon's avatar
Praetorius, Simon committed
209 210 211 212 213
				    PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
    }
214

Praetorius, Simon's avatar
Praetorius, Simon committed
215 216 217 218 219 220 221 222 223 224

    KSPSetType(kspSchur, KSPPREONLY);
    PC pcSub;
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell2);
    PCShellSetContext(pcSub, &matShellContext);

    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
225

Praetorius, Simon's avatar
Praetorius, Simon committed
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244

    // === Mass matrix solver ===

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    {
      Operator massOp(pressureFeSpace, pressureFeSpace);
      ZeroOrderTerm *massTerm = NULL;
      if ((!phase) || (*nu == 0.0))
	massTerm = new Simple_ZOT;
      else
        massTerm = new VecAtQP_ZOT(phase);
      massOp.addTerm(massTerm);
      massMatrix.assembleOperator(massOp);
      delete massTerm;
    }
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);

245

Praetorius, Simon's avatar
Praetorius, Simon committed
246 247 248
    if (neumannBC.size() > 0) {
      MSG("Number of neumann boundary indices: %d\n", neumannBC.size());
    }
249

Praetorius, Simon's avatar
Praetorius, Simon committed
250 251 252 253 254 255

    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    {
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
256
      SecondOrderTerm *laplaceTerm = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
257 258 259 260 261 262 263
      if ((!phase) || (*nu == 0.0))
	laplaceTerm = new Simple_SOT;
      else
	laplaceTerm = new VecAtQP_SOT(phase);
      laplaceOp.addTerm(laplaceTerm);
      for (size_t i = 0; i < neumannBC.size(); i++)
	addDirichletBC(laplaceMatrix, neumannBC[i]);
264

Praetorius, Simon's avatar
Praetorius, Simon committed
265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
      laplaceMatrix.assembleOperator(laplaceOp);
      delete laplaceTerm;
    }
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
    DOFVector<double> vz(pressureFeSpace, "vz");
    DOFVector<double> vp(pressureFeSpace, "vp");
    vx.interpol(solution->getDOFVector(0));
    if (dow >= 2)
      vy.interpol(solution->getDOFVector(1));
    if (dow >= 3)
      vz.interpol(solution->getDOFVector(2));

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);

    {
      Operator conDifOp(pressureFeSpace, pressureFeSpace);
288

Praetorius, Simon's avatar
Praetorius, Simon committed
289 290 291 292 293 294
      ZeroOrderTerm *conDif0 = NULL;
      SecondOrderTerm *conDif1 = NULL;
      FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;

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

Praetorius, Simon's avatar
Praetorius, Simon committed
296 297 298 299 300 301 302 303 304 305 306
	conDif0 = new Simple_ZOT(*invTau);
	conDifOp.addTerm(conDif0);
	conDif1 = new Simple_SOT(*nu);
	conDifOp.addTerm(conDif1);
	conDif2 = new VecAtQP_FOT(&vx, NULL, 0);
	conDifOp.addTerm(conDif2, GRD_PHI);
	if (dow >= 2) {
	  conDif3 = new VecAtQP_FOT(&vy, NULL, 1);
	  conDifOp.addTerm(conDif3, GRD_PHI);
	}
	if (dow == 3) {
307
	  conDif4 = new VecAtQP_FOT(&vz, NULL, 2);
Praetorius, Simon's avatar
Praetorius, Simon committed
308 309 310
	  conDifOp.addTerm(conDif4, GRD_PHI);
	}
      } else { // no phase given
311

Praetorius, Simon's avatar
Praetorius, Simon committed
312
	vp.interpol(phase);
313

Praetorius, Simon's avatar
Praetorius, Simon committed
314 315
	if (*nu > 0.0) {
	  conDif0 = new VecAtQP_ZOT(&vp, NULL, *invTau);
316
	  conDifOp.addTerm(conDif0);
Praetorius, Simon's avatar
Praetorius, Simon committed
317
	  conDif1 = new VecAtQP_SOT(&vp, NULL, *nu);
318
	  conDifOp.addTerm(conDif1);
Praetorius, Simon's avatar
Praetorius, Simon committed
319 320
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, NULL, 0);
	  conDifOp.addTerm(conDif2, GRD_PHI);
321

Praetorius, Simon's avatar
Praetorius, Simon committed
322 323 324 325 326 327 328 329
	  if (dow >= 2) {
	    conDif3 = new Vec2AtQP_FOT(&vy, &vp, NULL, 1);
	    conDifOp.addTerm(conDif3, GRD_PHI);
	  }
	  if (dow == 3) {
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, NULL, 2);
	    conDifOp.addTerm(conDif4, GRD_PHI);
	  }
330
	} else {
Praetorius, Simon's avatar
Praetorius, Simon committed
331
	  conDif0 = new VecAtQP_ZOT(&vp, new LinearInterpolation(*rho1,*rho2,*invTau));
332
	  conDifOp.addTerm(conDif0);
Praetorius, Simon's avatar
Praetorius, Simon committed
333
	  conDif1 = new VecAtQP_SOT(&vp, new LinearInterpolation(*nu1,*nu2));
334
	  conDifOp.addTerm(conDif1);
Praetorius, Simon's avatar
Praetorius, Simon committed
335 336
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, new LinearInterpolation2(*rho1,*rho2), 0);
	  conDifOp.addTerm(conDif2, GRD_PHI);
337

Praetorius, Simon's avatar
Praetorius, Simon committed
338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
	  if (dow >= 2) {
	    conDif3 = new Vec2AtQP_FOT(&vy, &vp, new LinearInterpolation2(*rho1,*rho2), 1);
	    conDifOp.addTerm(conDif3, GRD_PHI);
	  }
	  if (dow == 3) {
	    conDif4 = new Vec2AtQP_FOT(&vz, &vp, new LinearInterpolation2(*rho1,*rho2), 2);
	    conDifOp.addTerm(conDif4, GRD_PHI);
	  }
	}
      }/**/

      for (size_t i = 0; i < neumannBC.size(); i++)
	addDirichletBC(conDifMatrix, neumannBC[i]);
      conDifMatrix.assembleOperator(conDifOp);

      delete conDif0;
      delete conDif1;
      delete conDif2;
      if (dow >= 2)
	delete conDif3;
      if (dow == 3)
	delete conDif4;
    }
361

Praetorius, Simon's avatar
Praetorius, Simon committed
362 363 364 365 366 367 368 369 370

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


    // === Setup solver ===

    matShellContext.kspMass = massMatrixSolver->getSolver();
    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
371
    matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
Praetorius, Simon's avatar
Praetorius, Simon committed
372 373 374

    switch (massSolutionMode) {
    case 0:
375
      petsc_helper::setSolver(matShellContext.kspMass, "mass_",
Praetorius, Simon's avatar
Praetorius, Simon committed
376 377 378
			      KSPCG, PCJACOBI, 0.0, 1e-14, 2);
      break;
    case 1:
379
      petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON,
Praetorius, Simon's avatar
Praetorius, Simon committed
380 381 382 383 384 385 386 387
				    PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode);
    }

    switch (laplaceSolutionMode) {
    case 0:
388
      petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_",
Praetorius, Simon's avatar
Praetorius, Simon committed
389 390 391
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
      break;
    case 1:
392 393
      petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "laplace_",
				    KSPRICHARDSON, PCLU, MATSOLVERMUMPS,
Praetorius, Simon's avatar
Praetorius, Simon committed
394 395 396 397 398 399 400 401
				    0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode);
    }

    setConstantNullSpace(matShellContext.kspLaplace);

402
    MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n",
Praetorius, Simon's avatar
Praetorius, Simon committed
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
	t.elapsed());
  }


  void PetscSolverNavierStokes2::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes2::exitPreconditioner()");

    massMatrixSolver->destroyMatrixData();
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyMatrixData();
    laplaceMatrixSolver->destroyVectorData();
    conDifMatrixSolver->destroyMatrixData();
    conDifMatrixSolver->destroyVectorData();

    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    conDifMatrixSolver->destroyVectorData();
421 422


Praetorius, Simon's avatar
Praetorius, Simon committed
423 424 425 426 427 428 429 430 431 432
    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
} }