From 054249dcf1137d6e7d748726e21dc4479df11c92 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 23 Oct 2013 16:30:54 +0000 Subject: [PATCH] preconditioners added --- extensions/preconditioner/CahnHilliard_.cc | 217 +++++++++ extensions/preconditioner/CahnHilliard_.h | 66 +++ .../PetscSolverNavierStokes2.cc | 428 ++++++++++++++++++ .../preconditioner/PetscSolverNavierStokes2.h | 191 ++++++++ extensions/preconditioner/PetscSolverPfc.cc | 237 ++++++++++ extensions/preconditioner/PetscSolverPfc.h | 87 ++++ .../preconditioner/PhaseFieldCrystal_.cc | 87 ++++ .../preconditioner/PhaseFieldCrystal_.h | 78 ++++ 8 files changed, 1391 insertions(+) create mode 100644 extensions/preconditioner/CahnHilliard_.cc create mode 100644 extensions/preconditioner/CahnHilliard_.h create mode 100644 extensions/preconditioner/PetscSolverNavierStokes2.cc create mode 100644 extensions/preconditioner/PetscSolverNavierStokes2.h create mode 100644 extensions/preconditioner/PetscSolverPfc.cc create mode 100644 extensions/preconditioner/PetscSolverPfc.h create mode 100644 extensions/preconditioner/PhaseFieldCrystal_.cc create mode 100644 extensions/preconditioner/PhaseFieldCrystal_.h diff --git a/extensions/preconditioner/CahnHilliard_.cc b/extensions/preconditioner/CahnHilliard_.cc new file mode 100644 index 00000000..795b233c --- /dev/null +++ b/extensions/preconditioner/CahnHilliard_.cc @@ -0,0 +1,217 @@ +/****************************************************************************** + * + * 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 "CahnHilliard_.h" +// #include "Views.h" +#include "SignedDistFunctors.h" +#include "PhaseFieldConvert.h" + +#include "HL_SignedDistTraverse.h" +#include "Recovery.h" + +using namespace AMDiS; + +CahnHilliard_::CahnHilliard_(const std::string &name_) : + super(name_), + useMobility(false), + useReinit(false), + doubleWell(0), + gamma(1.0), + eps(0.1), + minusEps(-0.1), + epsInv(10.0), + minusEpsInv(-10.0), + epsSqr(0.01), + minusEpsSqr(-0.01) +{ + // parameters for CH + Parameters::get(name_ + "->use mobility", useMobility); // mobility + Parameters::get(name_ + "->gamma", gamma); // mobility + Parameters::get(name_ + "->epsilon", eps); // interface width + + // type of double well: 0= [0,1], 1= [-1,1] + Parameters::get(name_ + "->double-well type", doubleWell); + + Parameters::get(name + "->use reinit", useReinit); + + // transformation of the parameters + minusEps = -eps; + epsInv = 1.0/eps; + minusEpsInv = -epsInv; + epsSqr = sqr(eps); + minusEpsSqr = -epsSqr; +} + + +void CahnHilliard_::solveInitialProblem(AdaptInfo *adaptInfo) +{ + Flag initFlag = initDataFromFile(adaptInfo); + + if (!initFlag.isSet(DATA_ADOPTED)) { + int initialInterface = 0; + Initfile::get(name + "->initial interface", initialInterface); + double initialEps = eps; + Initfile::get(name + "->initial epsilon", initialEps); + + if (initialInterface == 0) { + /// horizontale Linie + double a= 0.0, dir= -1.0; + Initfile::get(name + "->line->pos", a); + Initfile::get(name + "->line->direction", dir); + prob->getSolution()->getDOFVector(1)->interpol(new Plane(a, dir)); + } + else if (initialInterface == 1) { + /// schraege Linie + double theta = m_pi/4.0; + prob->getSolution()->getDOFVector(1)->interpol(new PlaneRotation(0.0, theta, 1.0)); + transformDOFInterpolation(prob->getSolution()->getDOFVector(1),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min); + } + else if (initialInterface == 2) { + /// Ellipse + double a= 1.0, b= 1.0; + Initfile::get(name + "->ellipse->a", a); + Initfile::get(name + "->ellipse->b", b); + prob->getSolution()->getDOFVector(1)->interpol(new Ellipse(a,b)); + } + else if (initialInterface == 3) { + /// zwei horizontale Linien + double a= 0.75, b= 0.375; + Initfile::get(name + "->lines->pos1", a); + Initfile::get(name + "->lines->pos2", b); + prob->getSolution()->getDOFVector(1)->interpol(new Plane(a, -1.0)); + transformDOFInterpolation(prob->getSolution()->getDOFVector(1),new Plane(b, 1.0), new AMDiS::Max); + } + else if (initialInterface == 4) { + /// Kreis + double radius= 1.0; + Initfile::get(name + "->kreis->radius", radius); + prob->getSolution()->getDOFVector(1)->interpol(new Circle(radius)); + } else if (initialInterface == 5) { + /// Rechteck + double width = 0.5; + double height = 0.3; + WorldVector center; center.set(0.5); + Initfile::get(name + "->rectangle->width", width); + Initfile::get(name + "->rectangle->height", height); + Initfile::get(name + "->rectangle->center", center); + prob->getSolution()->getDOFVector(1)->interpol(new Rectangle(width, height, center)); + } + + // TODO: Redistancing einfügen! + if (useReinit) { + FiniteElemSpace* feSpace = FiniteElemSpace::provideFeSpace( + const_cast(getMesh()->getVertexAdmin()), + Lagrange::getLagrange(getMesh()->getDim(), 1), + getMesh(), + "P1"); + DOFVector tmp(feSpace, "tmp"); + tmp.interpol(prob->getSolution()->getDOFVector(1)); + + HL_SignedDistTraverse reinit("reinit", getMesh()->getDim()); + reinit.calcSignedDistFct(adaptInfo, &tmp); + +#ifndef HAVE_PARALLEL_DOMAIN_AMDIS + Recovery recovery(L2_NORM, 1); + recovery.recoveryUh(&tmp, *prob->getSolution()->getDOFVector(1)); +#else + prob->getSolution()->getDOFVector(1)->interpol(&tmp); +#endif + } + + + /// create phase-field from signed-dist-function + if (doubleWell == 0) { + forEachDOF(prob->getSolution()->getDOFVector(1), + new SignedDistToPhaseField(initialEps)); + } else { + forEachDOF(prob->getSolution()->getDOFVector(1), + new SignedDistToCh(initialEps)); + } + } +} + + +void CahnHilliard_::fillOperators() +{ + const FiniteElemSpace* feSpace = prob->getFeSpace(); + int degree = feSpace->getBasisFcts()->getDegree(); + + // c + Operator *opChMnew = new Operator(feSpace, feSpace); + opChMnew->addTerm(new Simple_ZOT); + Operator *opChMold = new Operator(feSpace, feSpace); + opChMold->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1))); + // -nabla*(grad(c)) + Operator *opChL = new Operator(feSpace, feSpace); + opChL->addTerm(new Simple_SOT); + + // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2 + Operator *opChLM = new Operator(feSpace, feSpace); + if (useMobility) { + if (doubleWell == 0) + opChLM->addTerm(new VecAtQP_SOT( + prob->getSolution()->getDOFVector(1), + new MobilityCH0(gamma, degree))); + else + opChLM->addTerm(new VecAtQP_SOT( + prob->getSolution()->getDOFVector(1), + new MobilityCH1(gamma, degree))); + } else + opChLM->addTerm(new Simple_SOT(gamma)); + + // -2*c_old^3 + 3/2*c_old^2 + Operator *opChMPowExpl = new Operator(feSpace, feSpace); + opChMPowExpl->addTerm(new VecAtQP_ZOT( + prob->getSolution()->getDOFVector(1), + new AMDiS::Pow<3>(-2.0, 3*degree))); + if (doubleWell == 0) { + opChMPowExpl->addTerm(new VecAtQP_ZOT( + prob->getSolution()->getDOFVector(1), + new AMDiS::Pow<2>(3.0/2.0, 2*degree))); + } + + // -3*c_old^2 * c + Operator *opChMPowImpl = new Operator(feSpace, feSpace); + opChMPowImpl->addTerm(new VecAtQP_ZOT( + prob->getSolution()->getDOFVector(1), + new AMDiS::Pow<2>(-3.0, 2*degree))); + if (doubleWell == 0) { + opChMPowImpl->addTerm(new VecAtQP_ZOT( + prob->getSolution()->getDOFVector(1), + NULL, 3.0)); + opChMPowImpl->addTerm(new Simple_ZOT(-0.5)); + } else { + opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0)); + } + + // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC] + // ---------------------------------------------------------------------- + prob->addMatrixOperator(*opChMPowImpl,0,1); /// < -3*c*c_old^2 , psi > + prob->addMatrixOperator(*opChL,0,1, &minusEpsSqr); /// < -eps^2*grad(c) , grad(psi) > + prob->addMatrixOperator(*opChMnew,0,0); /// < mu , psi > + // . . . vectorOperators . . . . . . . . . . . . . . . + prob->addVectorOperator(*opChMPowExpl,0); /// < -2*c_old^3 , psi > + + + // dt(c) = laplace(mu) - u*grad(c) + // ----------------------------------- + prob->addMatrixOperator(*opChMnew,1,1); /// < c , psi > + prob->addMatrixOperator(*opChLM,1,0, getTau()); /// < tau*grad(mu) , grad(psi) > + // . . . vectorOperators . . . . . . . . . . . . . . . + prob->addVectorOperator(*opChMold,1); /// < c^old , psi > + +} + diff --git a/extensions/preconditioner/CahnHilliard_.h b/extensions/preconditioner/CahnHilliard_.h new file mode 100644 index 00000000..9da54bc0 --- /dev/null +++ b/extensions/preconditioner/CahnHilliard_.h @@ -0,0 +1,66 @@ +/****************************************************************************** + * + * 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. + * + ******************************************************************************/ + +#ifndef CAHN_HILLIARD_PRECON_H +#define CAHN_HILLIARD_PRECON_H + +#include "AMDiS.h" +#include "BaseProblem.h" +#include "chns.h" + +using namespace AMDiS; + +class CahnHilliard_ : public BaseProblem +{ +public: // definition of types + + typedef BaseProblem super; + +public: // public methods + + CahnHilliard_(const std::string &name_); + ~CahnHilliard_() {}; + + void solveInitialProblem(AdaptInfo *adaptInfo); + + double getEpsilon() { return eps; } + int getDoubleWellType() { return doubleWell; } + + void fillOperators() override; + void fillBoundaryConditions() override {} + +protected: // protected variables + + bool useMobility; + bool useReinit; + + unsigned dim; + + int doubleWell; + + double gamma; + double eps; + double minusEps; + double epsInv; + double minusEpsInv; + double epsSqr; + double minusEpsSqr; +}; + + + +#endif // CAHN_HILLIARD_PRECON_H diff --git a/extensions/preconditioner/PetscSolverNavierStokes2.cc b/extensions/preconditioner/PetscSolverNavierStokes2.cc new file mode 100644 index 00000000..5792dd75 --- /dev/null +++ b/extensions/preconditioner/PetscSolverNavierStokes2.cc @@ -0,0 +1,428 @@ +/****************************************************************************** + * + * 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 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. + * + ******************************************************************************/ + + + +#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(ctx); + + // Project out constant null space + { + int vecSize; + VecGetSize(x, &vecSize); + PetscScalar vecSum; + VecSum(x, &vecSum); + vecSum = vecSum / static_cast(-1.0 * vecSize); + VecShift(x, vecSum); + } + + KSPSolve(data->kspLaplace, x, y); + MatMult(data->matConDif, y, x); + KSPSolve(data->kspMass, x, y); + + PetscFunctionReturn(0); + } + + + 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) + { + Parameters::get(name + "->navierstokes->pressure component", + 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"); + + Parameters::get(name + "->navierstokes->use old initial guess", + useOldInitialGuess); + + Parameters::get(name + "->navierstokes->velocity solver", + velocitySolutionMode); + + Parameters::get(name + "->navierstokes->mass solver", + massSolutionMode); + + Parameters::get(name + "->navierstokes->laplace solver", + laplaceSolutionMode); + + + Parameters::get(name + "->navierstokes->neumann boundary indices", + neumannBC); + } + + int PetscSolverNavierStokes2::solveLinearSystem(const SolverMatrix >& A, + SystemVector& x, + SystemVector& b, + bool createMatrixData, + bool storeMatrixData) + { + FUNCNAME("PetscSolverNavierStokes2::solveLinearSystem()"); + + 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); + + INFO(info, 8)("creation of parallel data structures needed %.5f seconds\n", + t.elapsed()); + + solvePetscMatrix(x, NULL); + + if (!storeMatrixData) { + destroyVectorData(); + destroyMatrixData(); + } + + return 0; + } + + + void PetscSolverNavierStokes2::initSolver(KSP &ksp) + { + // Create FGMRES based outer solver + KSPCreate(domainComm, &ksp); + KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); + 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, "ns_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations()); + + // Create null space information. + if (pressureNullSpace) + setConstantNullSpace(ksp, pressureComponent, true); + + if (useOldInitialGuess) + KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); + } + + void PetscSolverNavierStokes2::addDirichletBC(DOFMatrix& m, BoundaryType b) + { + DirichletBC *dirichletApply = + new DirichletBC(b, new AMDiS::Const >(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 velocityComponents; + for (size_t i = 0; i < static_cast(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) { + case 0: + petsc_helper::setSolver(kspVelocity, "", + KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); + break; + case 1: + petsc_helper::setSolverWithLu(kspVelocity, "", KSPPREONLY, + PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1); + break; + default: + ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode); + } + + + KSPSetType(kspSchur, KSPPREONLY); + PC pcSub; + KSPGetPC(kspSchur, &pcSub); + PCSetType(pcSub, PCSHELL); + PCShellSetApply(pcSub, pcSchurShell2); + PCShellSetContext(pcSub, &matShellContext); + + if (pressureNullSpace) + setConstantNullSpace(kspSchur); + + + // === 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); + + + if (neumannBC.size() > 0) { + MSG("Number of neumann boundary indices: %d\n", neumannBC.size()); + } + + + // === Laplace matrix solver === + + DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace); + { + Operator laplaceOp(pressureFeSpace, pressureFeSpace); + SecondOrderTerm *laplaceTerm = NULL; + 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]); + + laplaceMatrix.assembleOperator(laplaceOp); + delete laplaceTerm; + } + laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_")); + laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix); + + + // === Create convection-diffusion operator === + + DOFVector vx(pressureFeSpace, "vx"); + DOFVector vy(pressureFeSpace, "vy"); + DOFVector vz(pressureFeSpace, "vz"); + DOFVector 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); + + ZeroOrderTerm *conDif0 = NULL; + SecondOrderTerm *conDif1 = NULL; + FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL; + + if (!phase) { + MSG("INIT WITHOUT PHASE!\n"); + + 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) { + conDif4 = new VecAtQP_FOT(&vz, NULL, 2); + conDifOp.addTerm(conDif4, GRD_PHI); + } + } else { // no phase given + + vp.interpol(phase); + + if (*nu > 0.0) { + conDif0 = new VecAtQP_ZOT(&vp, NULL, *invTau); + conDifOp.addTerm(conDif0); + conDif1 = new VecAtQP_SOT(&vp, NULL, *nu); + conDifOp.addTerm(conDif1); + conDif2 = new Vec2AtQP_FOT(&vx, &vp, NULL, 0); + conDifOp.addTerm(conDif2, GRD_PHI); + + 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); + } + } else { + conDif0 = new VecAtQP_ZOT(&vp, new LinearInterpolation(*rho1,*rho2,*invTau)); + conDifOp.addTerm(conDif0); + conDif1 = new VecAtQP_SOT(&vp, new LinearInterpolation(*nu1,*nu2)); + conDifOp.addTerm(conDif1); + conDif2 = new Vec2AtQP_FOT(&vx, &vp, new LinearInterpolation2(*rho1,*rho2), 0); + conDifOp.addTerm(conDif2, GRD_PHI); + + 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; + } + + + conDifMatrixSolver = createSubSolver(pressureComponent, "condif_"); + conDifMatrixSolver->fillPetscMatrix(&conDifMatrix); + + + // === Setup solver === + + matShellContext.kspMass = massMatrixSolver->getSolver(); + matShellContext.kspLaplace = laplaceMatrixSolver->getSolver(); + matShellContext.matConDif = conDifMatrixSolver->getMatInterior(); + + switch (massSolutionMode) { + case 0: + petsc_helper::setSolver(matShellContext.kspMass, "mass_", + KSPCG, PCJACOBI, 0.0, 1e-14, 2); + break; + case 1: + petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON, + PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1); + break; + default: + ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode); + } + + switch (laplaceSolutionMode) { + case 0: + petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", + KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); + break; + case 1: + petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "laplace_", + KSPRICHARDSON, PCLU, MATSOLVERMUMPS, + 0.0, 1e-14, 1); + break; + default: + ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode); + } + + setConstantNullSpace(matShellContext.kspLaplace); + + MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n", + 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(); + + + delete massMatrixSolver; + massMatrixSolver = NULL; + + delete laplaceMatrixSolver; + laplaceMatrixSolver = NULL; + + delete conDifMatrixSolver; + conDifMatrixSolver = NULL; + } +} } diff --git a/extensions/preconditioner/PetscSolverNavierStokes2.h b/extensions/preconditioner/PetscSolverNavierStokes2.h new file mode 100644 index 00000000..091ff96d --- /dev/null +++ b/extensions/preconditioner/PetscSolverNavierStokes2.h @@ -0,0 +1,191 @@ +/****************************************************************************** + * + * 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 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. + * + ******************************************************************************/ + + + +/** \file PetscSolverNavierStokes2.h */ + +#ifndef AMDIS_PETSC_SOLVER_NAVIER_STOKES2_H +#define AMDIS_PETSC_SOLVER_NAVIER_STOKES2_H + +#include "parallel/PetscSolverGlobalMatrix.h" + +namespace AMDiS { namespace Parallel { + + using namespace std; + + struct NavierStokesSchurData2 { + KSP kspMass; + + KSP kspLaplace; + + Mat matConDif; + }; + + class PetscSolverNavierStokes2 : public PetscSolverGlobalMatrix + { + private: + + class LinearInterpolation : public AbstractFunction + { + public: + LinearInterpolation(double c1, double c2, double factor=1.0) + : AbstractFunction(0) + { a = (c1-c2)/2.0*factor; b = (c1+c2)/2.0*factor; cmin=std::min(c1,c2)*factor; cmax=std::max(c1,c2)*factor;} + + double operator()(const double& x) const + { double result = b+a*x; + if (resultcmax) result = cmax; + return result; + } + private: + double a,b,cmin,cmax; + }; + + + class LinearInterpolation2 : public BinaryAbstractFunction + { + public: + LinearInterpolation2(double c1, double c2, double factor=1.0) + : BinaryAbstractFunction(0) + { a = (c1-c2)/2.0*factor; b = (c1+c2)/2.0*factor; cmin=std::min(c1,c2)*factor; cmax=std::max(c1,c2)*factor;} + + double operator()(const double& u, const double& x) const + { + double result = b+a*x; + if (resultcmax) result = cmax; + return result * u; + } + private: + double a,b,cmin,cmax; + }; + + class EinsMinus : public AbstractFunction + { + public: + EinsMinus(double d) + : AbstractFunction(0), + c(d) + {} + + double operator()(const double& x) const + { + return c * std::max(1.0-x,0.000001); + } + private: + double c; + }; + + + + public: + /// Creator class + class Creator : public LinearSolverCreator + { + public: + virtual ~Creator() {} + + /// Returns a new PetscSolver object. + LinearSolver* create() + { + return new PetscSolverNavierStokes2(this->name); + } + }; + + PetscSolverNavierStokes2(string name); + + void setStokesData(double *invTauPtr, SystemVector *vec, double *nu1_=NULL, double *nu2_=NULL, double *rho1_=NULL, double *rho2_=NULL) + { + nu = new double; + (*nu) = 0.0; + invTau = invTauPtr; + solution = vec; + nu1=nu1_; + nu2=nu2_; + rho1=rho1_; + rho2=rho2_; + } + + void setStokesData(double *nuPtr, double *invTauPtr, SystemVector *vec) + { + nu = nuPtr; + invTau = invTauPtr; + solution = vec; + } + + + void setPhase(DOFVector *d) + { + phase = d; + } + + protected: + + /// Implementation of \ref LinearSolver::solveLinearSystem() + int solveLinearSystem(const SolverMatrix >& A, + SystemVector& x, + SystemVector& b, + bool createMatrixData, + bool storeMatrixData); + + void initSolver(KSP &ksp); + + void initPreconditioner(PC pc); + + void exitPreconditioner(PC pc); + + void addDirichletBC(DOFMatrix &m, BoundaryType b); + + private: + int pressureComponent; + + bool pressureNullSpace; + + /// If true, old solution is used for initial guess in solver phase. + bool useOldInitialGuess; + + /// 0: approximate solve 1: direct solver + int velocitySolutionMode; + + /// 0: approximate solve 1: direct solver + int massSolutionMode; + + /// 0: approximate solve 1: direct solver + int laplaceSolutionMode; + + PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver; + + NavierStokesSchurData2 matShellContext; + + double *nu, *invTau, *nu1,*nu2,*rho1,*rho2; + + SystemVector* solution; + + DOFVector* phase; + + std::vector neumannBC; + const Matrix* systemMat; + }; + +} } + +#endif diff --git a/extensions/preconditioner/PetscSolverPfc.cc b/extensions/preconditioner/PetscSolverPfc.cc new file mode 100644 index 00000000..93643707 --- /dev/null +++ b/extensions/preconditioner/PetscSolverPfc.cc @@ -0,0 +1,237 @@ +/****************************************************************************** + * + * 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.h" +#include "parallel/PetscHelper.h" +#include "parallel/PetscSolverGlobalMatrix.h" + +namespace AMDiS { namespace Parallel { + + using namespace std; + PetscErrorCode pcPfcMatMultSchur(Mat mat, Vec b, Vec x) + { + void *ctx; + MatShellGetContext(mat, &ctx); + PfcData* data = static_cast(ctx); + + Vec y1, y2, y3; + VecDuplicate(b, &y1); + VecDuplicate(b, &y2); + VecDuplicate(b, &y3); + + MatMult(data->matK, b, y1); // y1 := K*b + KSPSolve(data->kspM, y1, y2); // M*y2 = y1 + MatMult(data->matK, y2, y3); // y3 := K*y2 +// KSPSolve(data->kspM, y3, y2); // M*y2 = y3 +// MatMult(data->matK, y2, y3); // y3 := K*y2 + + VecScale(y3, data->delta); // y3 = delta*y3 + MatMultAdd(data->matM, b, y3, x); // x = M*b + y3 + VecAXPY(x, -2.0*(data->delta), y1); // x := M*b - 2*delta*K*b + delta*B*b + + VecDestroy(&y1); + VecDestroy(&y2); + VecDestroy(&y3); + + PetscFunctionReturn(0); + } + + /// solve Pfc Preconditioner + PetscErrorCode pcPfcShell(PC pc, Vec b, Vec x) // solve Px=b + { FUNCNAME("pcPfcShell()"); + + void *ctx; + PCShellGetContext(pc, &ctx); + PfcData* 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->kspSchur, tmp, x2); // S*x2 = tmp + 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->kspM2, tmp, x3); // M*x3 = tmp + + VecDestroy(&y1); + VecDestroy(&y2); + VecDestroy(&tmp); + + PetscFunctionReturn(0); + } + + void PetscSolverPfc::initSolver(KSP &ksp) + { + // Create FGMRES based outer solver + KSPCreate(meshDistributor->getMpiComm(0), &ksp); + KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); + 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 PetscSolverPfc::initPreconditioner(PC pc) + { + FUNCNAME("PetscSolverPfc::initPreconditioner()"); + + TEST_EXIT(tau)("tau pointer not set!\n"); + + PCSetType(getPc(), PCSHELL); + PCShellSetApply(getPc(), pcPfcShell); + PCShellSetContext(getPc(), &data); + + 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(); + + solverM2 = createSubSolver(0, "M2_"); + solverM2->fillPetscMatrix(&matrixM); + data.kspM2 = solverM2->getSolver(); + + // create laplace-matrix + DOFMatrix matrixMpK(feSpace, feSpace); + Operator laplaceOp2(feSpace, feSpace); + Simple_SOT sot2(sqrt(*tau)); + laplaceOp2.addTerm(&zot); + laplaceOp2.addTerm(&sot2); + matrixMpK.assembleOperator(laplaceOp2); + + solverMpK = createSubSolver(0, "MpK_"); + solverMpK->fillPetscMatrix(&matrixMpK); + matMpK = solverMpK->getMatInterior(); + data.kspMpK = solverMpK->getSolver(); + + // 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(); + + // === Setup preconditioner data === + data.delta = sqrt(*tau); + data.tau = (*tau); + + petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCJACOBI, 0.0, 1e-14, 2); + petsc_helper::setSolver(data.kspM2, "M2_", KSPCG, PCBJACOBI, 0.0, 1e-14, 5); + petsc_helper::setSolver(data.kspMpK, "MpK_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); + + PetscInt mLocal, mGlobal; + MatGetLocalSize(data.matM, &mLocal, PETSC_NULL); + MatGetSize(data.matM, &mGlobal, PETSC_NULL); + MatCreateShell(meshDistributor->getMpiComm(0), + mLocal, mLocal, + mGlobal, mGlobal, + &data, + &matSchur); + MatShellSetOperation(matSchur, MATOP_MULT, (void(*)(void))pcPfcMatMultSchur); + + KSPCreate(meshDistributor->getMpiComm(0), &data.kspSchur); + KSPSetOperators(data.kspSchur, matSchur, matMpK, SAME_NONZERO_PATTERN); + petsc_helper::setSolver(data.kspSchur, "schur_", KSPFGMRES, PCJACOBI, 1e-6, 1e-8, 5); + KSPSetFromOptions(data.kspSchur); + } + + + void PetscSolverPfc::exitPreconditioner(PC pc) + { + FUNCNAME("PetscSolverPfc::exit()"); + + MatDestroy(&matMpK); + MatDestroy(&matSchur); + + solverM->destroyMatrixData(); + solverK->destroyMatrixData(); + solverMpK->destroyMatrixData(); + + solverM->destroyVectorData(); + solverK->destroyVectorData(); + solverMpK->destroyVectorData(); + + delete solverM; + solverM = NULL; + delete solverK; + solverK = NULL; + delete solverMpK; + solverMpK = NULL; + } + + + PetscSolver* PetscSolverPfc::createSubSolver(int component, + string kspPrefix) + { + FUNCNAME("PetscSolverPfc::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; + } +} } diff --git a/extensions/preconditioner/PetscSolverPfc.h b/extensions/preconditioner/PetscSolverPfc.h new file mode 100644 index 00000000..80d0477d --- /dev/null +++ b/extensions/preconditioner/PetscSolverPfc.h @@ -0,0 +1,87 @@ +/****************************************************************************** + * + * 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. + * + ******************************************************************************/ + +#ifndef AMDIS_PETSC_PRECON_PFC_H +#define AMDIS_PETSC_PRECON_PFC_H + +#include "parallel/PetscSolverGlobalBlockMatrix.h" + +namespace AMDiS { namespace Parallel { + + struct PfcData { + KSP kspM, kspM2, kspMpK, kspSchur; + Mat matM, matK; + double delta, tau; + }; + + class PetscSolverPfc : public PetscSolverGlobalBlockMatrix + { + public: + /// Creator class + class Creator : public LinearSolverCreator + { + public: + virtual ~Creator() {} + + /// Returns a new PetscSolver object. + LinearSolver* create() + { + return new PetscSolverPfc(this->name); + } + }; + + PetscSolverPfc(std::string name) + : PetscSolverGlobalBlockMatrix(name), + solverM(NULL), + solverM2(NULL), + solverK(NULL), + solverMpK(NULL), + useOldInitialGuess(false), + tau(NULL) + { + Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess); + } + + void setData(double *tauPtr) + { + tau = tauPtr; + } + + protected: + void initSolver(KSP &ksp); + + void initPreconditioner(PC pc); + void exitPreconditioner(PC pc); + + PetscSolver* createSubSolver(int component, string kspPrefix); + + private: + Mat matMpK; + Mat matSchur; + PfcData data; + + PetscSolver *solverM, *solverM2, *solverK, *solverMpK; + + bool useOldInitialGuess; + SystemVector* solution; + + double *tau; + }; + +} } + +#endif // AMDIS_PETSC_PRECON_PFC_H diff --git a/extensions/preconditioner/PhaseFieldCrystal_.cc b/extensions/preconditioner/PhaseFieldCrystal_.cc new file mode 100644 index 00000000..5185c005 --- /dev/null +++ b/extensions/preconditioner/PhaseFieldCrystal_.cc @@ -0,0 +1,87 @@ +#include "PhaseFieldCrystal_.h" + +using namespace std; +using namespace AMDiS; + +PhaseFieldCrystal_::PhaseFieldCrystal_(const std::string &name_, bool createProblem) : + super(name_, createProblem), + useMobility(false), + tempParameter(-0.6), + r(-0.4), // temperature deviation + rho0(1.0), // liquid density + density(-0.3), // mean density + two(2.0), + minus2(-2.0) +{ + Parameters::get(name + "->r",r); + Parameters::get(name + "->rho0", rho0); + Parameters::get(name + "->density", density); + Parameters::get(name + "->use mobility", useMobility); + tempParameter= -(1.0+r); +} + + +void PhaseFieldCrystal_::fillOperators() +{ + + int degree = prob->getFeSpace(0)->getBasisFcts()->getDegree(); + + // phi*rho + Operator *opMnew = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + opMnew->addTerm(new Simple_ZOT); + Operator *opMold = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + opMold->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1), NULL)); + + // -nabla*(phi*grad(rho)) + Operator *opL = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + opL->addTerm(new Simple_SOT); + + Operator *opLM = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + double M0 = 1.0; + Parameters::get(name + "->M0", M0); + if (useMobility) // non-constant mobility + opLM->addTerm(new VecAtQP_SOT( + prob->getSolution()->getDOFVector(1), + new MobilityPfc(density, M0, degree))); + else // constant mobility + opLM->addTerm(new Simple_SOT(M0)); + + // -(1+r)*phi*rho + Operator *opMTemp = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + opMTemp->addZeroOrderTerm(new Simple_ZOT()); + // -2*rho_old^3 + Operator *opMPowExpl = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + //opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new Pow3Functor(-2.0))); + opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new AMDiS::Pow<3>(1.0, 3*degree))); + // -3*rho_old^2 + Operator *opMPowImpl = new Operator(prob->getFeSpace(0), prob->getFeSpace(0)); + opMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new AMDiS::Pow<2>(-3.0, 2*degree))); + + + // dt(rho) = laplace(mu) - u*grad(rho) + // ----------------------------------- + prob->addMatrixOperator(opMnew, 1, 1); + prob->addMatrixOperator(opLM, 1, 0, getTau(), getTau()); // -laplace(mu) + // . . . vectorOperators . . . . . . . . . . . . . . . + prob->addVectorOperator(opMold, 1); + + // mu-2*nu-laplace(nu)-(1+r)*rho = (rho_old^3) + ExtPot - eps^2/(rho_old+0.9) + // ---------------------------------------------------------------------- + prob->addMatrixOperator(opMTemp, 0, 1, getTempParameter(), getTempParameter()); // -phi*(1+r)*rho + //prob->addMatrixOperator(opMPowImpl, 0, 1); // -3*rho*rho_old^2 + prob->addMatrixOperator(opL, 0, 1, &two, &two); // -2*phi*laplace(rho) * psi + prob->addMatrixOperator(opMnew, 0, 0); // phi*mu * psi + prob->addMatrixOperator(opL, 0, 2); // phi*grad(nu) * grad(psi) + // prob.addMatrixOperator(opMnew,2,2,&minus2); + // . . . vectorOperators . . . . . . . . . . . . . . . + prob->addVectorOperator(opMPowExpl, 0); // -2*phi^old*rho_old^3 + + // 0 = nu-laplace(rho) + // ------------------- + prob->addMatrixOperator(opL, 2, 1); // -laplace(rho) + prob->addMatrixOperator(opMnew, 2, 2); // nu +} + + +void PhaseFieldCrystal_::fillBoundaryConditions() +{} diff --git a/extensions/preconditioner/PhaseFieldCrystal_.h b/extensions/preconditioner/PhaseFieldCrystal_.h new file mode 100644 index 00000000..5178f0ec --- /dev/null +++ b/extensions/preconditioner/PhaseFieldCrystal_.h @@ -0,0 +1,78 @@ +/** \file PhaseFieldCrystal_.h */ + +#ifndef PHASE_FIELD_CRYSTAL_PRECON_H_ +#define PHASE_FIELD_CRYSTAL_PRECON_H_ + +#include "AMDiS.h" +#include "BaseProblem.h" +#include "ExtendedProblemStat.h" + +using namespace AMDiS; + +/** Phase-field Crystal problem + */ +class PhaseFieldCrystal_ : public BaseProblem +{ +public: // typedefs + + typedef BaseProblem super; + +public: + + PhaseFieldCrystal_(const std::string &name_, bool createProblem = true); + ~PhaseFieldCrystal_() {} + + double *getTempParameter() { return &tempParameter; } + + virtual void fillOperators(); + virtual void fillBoundaryConditions(); + +protected: + + bool useMobility; + + double tempParameter; + double r; + double rho0; + double density; + double two; + double minus2; +}; + + +/** \ingroup MainInstat + * \brief + * Abstract function to calculate the pure PFC-Energy + */ +class Energy : public BinaryAbstractFunction +{ + public: + Energy() : BinaryAbstractFunction(4) { } + + double operator()(const double &rho, const double &mu) const { + return -0.25*sqr(sqr(rho)) + 0.5*rho*mu; } +}; + + +class MobilityPfc : public AbstractFunction +{ + public: + MobilityPfc(double density_ = -0.3, double factor_ = 1.0, int degree=1) : + AbstractFunction(degree+1), + density(density_), + factor(factor_), + delta(1.e-6) { } + + double operator()(const double &rho) const + { + double mobility= abs(rho + 1.5)*factor; + return std::max(mobility, 0.0); + } + + protected: + double density; + double factor; + double delta; +}; + +#endif // PHASE_FIELD_CRYSTAL_PRECON_H_ -- GitLab