PetscSolverNavierStokes2.cc 12.2 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
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
Praetorius, Simon's avatar
Praetorius, Simon committed
137
    if (getInfo() >= 10)
138
      petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
Praetorius, Simon's avatar
Praetorius, Simon committed
139
    else if (getInfo() >= 20)
140
      petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
Praetorius, Simon's avatar
Praetorius, Simon committed
141
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
142

Praetorius, Simon's avatar
Praetorius, Simon committed
143 144 145
    // Create null space information.
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
146

Praetorius, Simon's avatar
Praetorius, Simon committed
147 148 149
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }
150

Praetorius, Simon's avatar
Praetorius, Simon committed
151 152
  void PetscSolverNavierStokes2::addDirichletBC(DOFMatrix& m, BoundaryType b)
  {
153
    DirichletBC *dirichletApply =
Praetorius, Simon's avatar
Praetorius, Simon committed
154 155 156 157 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
      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) {
199 200
    case 0:
      petsc_helper::setSolver(kspVelocity, "",
Praetorius, Simon's avatar
Praetorius, Simon committed
201 202 203
 			      KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
      break;
    case 1:
204
      petsc_helper::setSolverWithLu(kspVelocity, "", KSPPREONLY,
Praetorius, Simon's avatar
Praetorius, Simon committed
205 206 207 208 209
				    PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
      break;
    default:
      ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
    }
210

Praetorius, Simon's avatar
Praetorius, Simon committed
211 212 213 214 215 216 217 218 219 220

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

    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
221

Praetorius, Simon's avatar
Praetorius, Simon committed
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240

    // === 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);

241

Praetorius, Simon's avatar
Praetorius, Simon committed
242 243 244
    if (neumannBC.size() > 0) {
      MSG("Number of neumann boundary indices: %d\n", neumannBC.size());
    }
245

Praetorius, Simon's avatar
Praetorius, Simon committed
246 247 248 249 250 251

    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    {
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
252
      SecondOrderTerm *laplaceTerm = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
253 254 255 256 257 258 259
      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]);
260

Praetorius, Simon's avatar
Praetorius, Simon committed
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283
      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);
284

Praetorius, Simon's avatar
Praetorius, Simon committed
285 286 287 288 289 290
      ZeroOrderTerm *conDif0 = NULL;
      SecondOrderTerm *conDif1 = NULL;
      FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;

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

Praetorius, Simon's avatar
Praetorius, Simon committed
292 293 294 295 296 297 298 299 300 301 302
	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) {
303
	  conDif4 = new VecAtQP_FOT(&vz, NULL, 2);
Praetorius, Simon's avatar
Praetorius, Simon committed
304 305 306
	  conDifOp.addTerm(conDif4, GRD_PHI);
	}
      } else { // no phase given
307

Praetorius, Simon's avatar
Praetorius, Simon committed
308
	vp.interpol(phase);
309

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

Praetorius, Simon's avatar
Praetorius, Simon committed
318 319 320 321 322 323 324 325
	  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);
	  }
326
	} else {
Praetorius, Simon's avatar
Praetorius, Simon committed
327
	  conDif0 = new VecAtQP_ZOT(&vp, new LinearInterpolation(*rho1,*rho2,*invTau));
328
	  conDifOp.addTerm(conDif0);
Praetorius, Simon's avatar
Praetorius, Simon committed
329
	  conDif1 = new VecAtQP_SOT(&vp, new LinearInterpolation(*nu1,*nu2));
330
	  conDifOp.addTerm(conDif1);
Praetorius, Simon's avatar
Praetorius, Simon committed
331 332
	  conDif2 = new Vec2AtQP_FOT(&vx, &vp, new LinearInterpolation2(*rho1,*rho2), 0);
	  conDifOp.addTerm(conDif2, GRD_PHI);
333

Praetorius, Simon's avatar
Praetorius, Simon committed
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
	  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;
    }
357

Praetorius, Simon's avatar
Praetorius, Simon committed
358 359 360 361 362 363 364 365 366

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


    // === Setup solver ===

    matShellContext.kspMass = massMatrixSolver->getSolver();
    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
367
    matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
Praetorius, Simon's avatar
Praetorius, Simon committed
368 369 370

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

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

    setConstantNullSpace(matShellContext.kspLaplace);

398
    MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n",
Praetorius, Simon's avatar
Praetorius, Simon committed
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416
	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();
417 418


Praetorius, Simon's avatar
Praetorius, Simon committed
419 420 421 422 423 424 425 426 427 428
    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;

    delete conDifMatrixSolver;
    conDifMatrixSolver = NULL;
  }
} }