PetscSolverNavierStokes.cc 5.01 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



#include "parallel/PetscSolverNavierStokes.h"
15
#include "parallel/PetscHelper.h"
16 17 18 19 20

namespace AMDiS {

  using namespace std;

21

22
  PetscErrorCode pcSchurShell(PC pc, Vec x, Vec y)
Thomas Witkowski's avatar
Thomas Witkowski committed
23 24 25
  {
    void *ctx;
    PCShellGetContext(pc, &ctx);
26
    NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
Thomas Witkowski's avatar
Thomas Witkowski committed
27

28 29 30
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
31
  }
32
  
33 34 35 36 37

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
      massMatrixSolver(NULL),
38 39 40 41
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
      solution(NULL)
42 43 44 45 46 47 48 49 50 51 52
  {
    Parameters::get(initFileStr + "->stokes->pressure component", 
		    pressureComponent);

    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");

    TEST_EXIT(pressureComponent == 2)("Fixed for pressure component = 2\n");
  }


53
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
54 55 56
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

57
    // Create FGMRES based outer solver
58
    KSPCreate(mpiCommGlobal, &ksp);
59
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
60
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
61
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
62

63 64
    // Create null space information.
    setConstantNullSpace(ksp, pressureComponent);
65 66 67 68 69 70 71
  }


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

72 73 74
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
75

76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);

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

    petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);

    KSPSetType(kspSchur, KSPPREONLY);
98
    PC pcSub;
99 100 101 102
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
103

104 105 106 107
    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
    KSPSetNullSpace(kspSchur, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
108

109 110

    // === Mass matrix solver ===
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
    Simple_ZOT zot;
    massOp.addTerm(&zot);
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
    Simple_SOT sot;
    laplaceOp.addTerm(&sot);
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

135 136 137 138
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
139 140 141

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
142
    Simple_ZOT conDif0(*invTau);
143
    conDifOp.addTerm(&conDif0);
144
    Simple_SOT conDif1(*nu);
145
    conDifOp.addTerm(&conDif1);
146
    VecAtQP_FOT conDif2(&vx, &idFct, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
147
    conDifOp.addTerm(&conDif2, GRD_PHI);
148
    VecAtQP_FOT conDif3(&vy, &idFct, 1);
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    conDifOp.addTerm(&conDif3, GRD_PHI);
150 151 152 153 154 155

    conDifMatrix.assembleOperator(conDifOp);

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

156 157

    // === Setup solver ===
158 159 160 161 162

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

163 164
    petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
    petsc_helper::setSolver(matShellContext.kspLaplace, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
165 166
  }

167

168
}