/****************************************************************************** * * Extension of AMDiS - Adaptive multidimensional simulations * * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis * * Authors: 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. * * * See also license.opensource.txt in the distribution. * ******************************************************************************/ #include "PetscSolverPfc_diag.h" #include "parallel/PetscHelper.h" #include "parallel/PetscSolverGlobalMatrix.h" namespace AMDiS { namespace Parallel { using namespace std; /// solve Pfc Preconditioner PetscErrorCode pcPfcDiagShell(PC pc, Vec b, Vec x) // solve Px=b { FUNCNAME("pcPfcDiagShell()"); void *ctx; PCShellGetContext(pc, &ctx); PfcDiagData* data = static_cast(ctx); Vec b1, b2, b3, x1, x2, x3; VecNestGetSubVec(b, 0, &b1); VecNestGetSubVec(b, 1, &b2); VecNestGetSubVec(b, 2, &b3); VecNestGetSubVec(x, 0, &x1); VecNestGetSubVec(x, 1, &x2); VecNestGetSubVec(x, 2, &x3); // Hilfsvariablen Vec y1, y2, tmp; VecDuplicate(b1, &y1); VecDuplicate(b1, &y2); VecDuplicate(b1, &tmp); KSPSolve(data->kspM, b1, y1); // M*y1 = b1 MatMult(data->matK, y1, tmp); // tmp := K*y1 VecAYPX(tmp, -(data->tau), b2); // tmp := b2 - tau*tmp KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp MatMult(data->matM, y2, tmp); // tmp := M*y2 KSPSolve(data->kspS, tmp, x2); VecCopy(x2, x1); // x1 := x2 VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1 MatMult(data->matK, x2, tmp); // tmp := K*x2 VecAYPX(tmp, -1.0, b3); // tmp := b3 - tmp KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp VecDestroy(&y1); VecDestroy(&y2); VecDestroy(&tmp); return 0; } void PetscSolverPfcDiag::initSolver(KSP &ksp) { // Create FGMRES based outer solver KSPCreate(meshDistributor->getMpiComm(0), &ksp); #if (PETSC_VERSION_MINOR >= 5) KSPSetOperators(ksp, getMatInterior(), getMatInterior()); #else KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); #endif if (getInfo() >= 10) KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL); else if (getInfo() >= 20) KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL); petsc_helper::setSolver(ksp, "pfc_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations()); KSPSetFromOptions(ksp); if (useOldInitialGuess) KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); } void PetscSolverPfcDiag::initPreconditioner(PC pc) { FUNCNAME("PetscSolverPfcDiag::initPreconditioner()"); TEST_EXIT(tau)("tau pointer not set!\n"); PCSetType(getPc(), PCSHELL); PCShellSetApply(getPc(), pcPfcDiagShell); PCShellSetContext(getPc(), &data); double delta = sqrt(*tau); const FiniteElemSpace *feSpace = componentSpaces[0]; // create mass-matrix DOFMatrix matrixM(feSpace, feSpace); Operator massOp(feSpace, feSpace); Simple_ZOT zot; massOp.addTerm(&zot); matrixM.assembleOperator(massOp); solverM = createSubSolver(0, "M_"); solverM->fillPetscMatrix(&matrixM); data.matM = solverM->getMatInterior(); data.kspM = solverM->getSolver(); // create MpK-matrix DOFMatrix matrixMpK(feSpace, feSpace); Operator laplaceOp2(feSpace, feSpace); Simple_SOT sot2(delta); laplaceOp2.addTerm(&zot); laplaceOp2.addTerm(&sot2); matrixMpK.assembleOperator(laplaceOp2); solverMpK = createSubSolver(0, "MpK_"); solverMpK->fillPetscMatrix(&matrixMpK); matMpK = solverMpK->getMatInterior(); data.kspMpK = solverMpK->getSolver(); // create MpK2-matrix // create laplace-matrix DOFMatrix matrixK(feSpace, feSpace); Operator laplaceOp(feSpace, feSpace); Simple_SOT sot; laplaceOp.addTerm(&sot); matrixK.assembleOperator(laplaceOp); solverK = createSubSolver(0, "K_"); solverK->fillPetscMatrix(&matrixK); data.matK = solverK->getMatInterior(); // create Matrix S solverS = createSubSolver(0, "MpK2_"); solverS->fillPetscMatrix(&matrixM); // just to create a dummy matrix matS = solverS->getMatInterior(); data.kspS = solverS->getSolver(); Vec x; // need to initialize vector solverS->createVec(*solverS->getDofMapping(), x); Mat DK; MatDuplicate(data.matK, MAT_COPY_VALUES, &DK); MatGetDiagonal(data.matM, x); VecReciprocal(x); MatDiagonalScale(DK, x, PETSC_NULL); // DK := M_D^(-1)*K MatMatMult(data.matK, DK, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matS); // S := K*DK MatAYPX(matS, delta, data.matM, DIFFERENT_NONZERO_PATTERN); // S = delta*S + M MatAXPY(matS, (-2.0*delta), data.matK, DIFFERENT_NONZERO_PATTERN); // S = S - 2*sqrt(delta)*K VecDestroy(&x); MatDestroy(&DK); #if (PETSC_VERSION_MINOR >= 5) KSPSetOperators(data.kspS, matS, matS); #else KSPSetOperators(data.kspS, matS, matS, SAME_NONZERO_PATTERN); #endif // === Setup preconditioner data === data.delta = delta; data.tau = (*tau); int nIterM=5, nIterMpL=20, nIterMpL2=10; double tolM=1.e-6, tolMpL=1.e-6, tolMpL2=1.e-6; Parameters::get("precon_pfc_M->max iteration", nIterM); Parameters::get("precon_pfc_MpL->max iteration", nIterMpL); Parameters::get("precon_pfc_MpL2->max iteration", nIterMpL2); Parameters::get("precon_pfc_M->tolerance", tolM); Parameters::get("precon_pfc_MpL->tolerance", tolMpL); Parameters::get("precon_pfc_MpL2->tolerance", tolMpL2); petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCBJACOBI, 0.0, tolM, nIterM); petsc_helper::setSolver(data.kspMpK, "MpK_", KSPCG, PCBJACOBI, 0.0, tolMpL, nIterMpL); petsc_helper::setSolver(data.kspS, "MpK2_", KSPCG, PCBJACOBI, 0.0, tolMpL2, nIterMpL2); } void PetscSolverPfcDiag::exitPreconditioner(PC pc) { FUNCNAME("PetscSolverPfcDiag::exit()"); // MatDestroy(&matMpK); // MatDestroy(&matMpK2); solverM->destroyMatrixData(); solverK->destroyMatrixData(); solverMpK->destroyMatrixData(); solverS->destroyMatrixData(); solverM->destroyVectorData(); solverK->destroyVectorData(); solverMpK->destroyVectorData(); solverS->destroyVectorData(); delete solverM; solverM = NULL; delete solverK; solverK = NULL; delete solverMpK; solverMpK = NULL; delete solverS; solverS = NULL; } PetscSolver* PetscSolverPfcDiag::createSubSolver(int component, string kspPrefix) { FUNCNAME("PetscSolverPfcDiag::createSubSolver()"); vector fe; fe.push_back(componentSpaces[component]); PetscSolver* subSolver = new PetscSolverGlobalMatrix(""); subSolver->setKspPrefix(kspPrefix.c_str()); subSolver->setMeshDistributor(meshDistributor, 0); subSolver->init(fe, fe); ParallelDofMapping &subDofMap = subSolver->getDofMap(); subDofMap[0] = dofMap[component]; subDofMap.update(); return subSolver; } } }