// // 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" #include "parallel/PetscHelper.h" #include "TransformDOF.h" namespace AMDiS { using namespace std; PetscErrorCode pcSchurShell(PC pc, Vec x, Vec y) { void *ctx; PCShellGetContext(pc, &ctx); NavierStokesSchurData* data = static_cast(ctx); KSPSolve(data->kspLaplace, x, y); MatMult(data->matConDif, y, x); KSPSolve(data->kspMass, x, y); } PetscSolverNavierStokes::PetscSolverNavierStokes(string name) : PetscSolverGlobalMatrix(name), pressureComponent(-1), massMatrixSolver(NULL), laplaceMatrixSolver(NULL), nu(NULL), invTau(NULL), solution(NULL), phase(NULL) { Parameters::get(initFileStr + "->navierstokes->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"); } void PetscSolverNavierStokes::initSolver(KSP &ksp) { FUNCNAME("PetscSolverNavierStokes::initSolver()"); // Create FGMRES based outer solver KSPCreate(mpiCommGlobal, &ksp); KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL); petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100); // Create null space information. setConstantNullSpace(ksp, pressureComponent, true); } void PetscSolverNavierStokes::initPreconditioner(PC pc) { FUNCNAME("PetscSolverNavierStokes::initPreconditioner()"); TEST_EXIT(nu)("nu pointer not set!\n"); TEST_EXIT(invTau)("invtau pointer not set!\n"); TEST_EXIT(solution)("solution pointer not set!\n"); vector velocityComponents; velocityComponents.push_back(0); velocityComponents.push_back(1); 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, "velocity_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); KSPSetType(kspSchur, KSPPREONLY); PC pcSub; KSPGetPC(kspSchur, &pcSub); PCSetType(pcSub, PCSHELL); PCShellSetApply(pcSub, pcSchurShell); PCShellSetContext(pcSub, &matShellContext); setConstantNullSpace(kspSchur); // === Mass matrix solver === const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent]; DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace); Operator massOp(pressureFeSpace, pressureFeSpace); // if (!phase) massOp.addTerm(new Simple_ZOT); // else // massOp.addTerm(new VecAtQP_ZOT(phase, &idFct)); massMatrix.assembleOperator(massOp); massMatrixSolver = createSubSolver(pressureComponent, "mass_"); massMatrixSolver->fillPetscMatrix(&massMatrix); // VtkWriter::writeFile(phase, "phase.vtu"); // === Laplace matrix solver === DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace); Operator laplaceOp(pressureFeSpace, pressureFeSpace); if (!phase) laplaceOp.addTerm(new Simple_SOT); else laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct)); laplaceMatrix.assembleOperator(laplaceOp); laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_"); laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix); // === Create convection-diffusion operator === DOFVector vx(pressureFeSpace, "vx"); DOFVector vy(pressureFeSpace, "vy"); DOFVector vp(pressureFeSpace, "vp"); vx.interpol(solution->getDOFVector(0)); vy.interpol(solution->getDOFVector(1)); DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace); Operator conDifOp(pressureFeSpace, pressureFeSpace); if (!phase) { MSG("INIT WITHOUT PHASE!\n"); Simple_ZOT *conDif0 = new Simple_ZOT(*invTau); conDifOp.addTerm(conDif0); Simple_SOT *conDif1 = new Simple_SOT(*nu); conDifOp.addTerm(conDif1); VecAtQP_FOT *conDif2 = new VecAtQP_FOT(&vx, &idFct, 0); conDifOp.addTerm(conDif2, GRD_PHI); VecAtQP_FOT *conDif3 = new VecAtQP_FOT(&vy, &idFct, 1); conDifOp.addTerm(conDif3, GRD_PHI); } else { MSG("INIT WITH PHASE!\n"); vp.interpol(phase); VecAtQP_ZOT *conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau)); conDifOp.addTerm(conDif0); VecAtQP_SOT *conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu)); conDifOp.addTerm(conDif1); Vec2AtQP_FOT *conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0); conDifOp.addTerm(conDif2, GRD_PHI); Vec2AtQP_FOT *conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1); conDifOp.addTerm(conDif3, GRD_PHI); } conDifMatrix.assembleOperator(conDifOp); conDifMatrixSolver = createSubSolver(pressureComponent, "condif_"); conDifMatrixSolver->fillPetscMatrix(&conDifMatrix); // === Setup solver === matShellContext.kspMass = massMatrixSolver->getSolver(); matShellContext.kspLaplace = laplaceMatrixSolver->getSolver(); matShellContext.matConDif = conDifMatrixSolver->getMatInterior(); petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2); // petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); KSPSetType(matShellContext.kspLaplace, KSPRICHARDSON); KSPSetTolerances(matShellContext.kspLaplace, 0.0, 1e-14, PETSC_DEFAULT, 1); KSPSetOptionsPrefix(matShellContext.kspLaplace, "laplace_"); KSPSetFromOptions(matShellContext.kspLaplace); MatNullSpace matNullSpace; MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace); KSPSetNullSpace(matShellContext.kspLaplace, matNullSpace); MatNullSpaceDestroy(&matNullSpace); { PC pc; KSPGetPC(matShellContext.kspLaplace, &pc); PCSetType(pc, PCHYPRE); PCSetFromOptions(pc); } // setConstantNullSpace(matShellContext.kspLaplace); } }