From 014e5191679f22f06168402353b3e857328e6cf5 Mon Sep 17 00:00:00 2001
From: Sebastian Aland <sebastian.aland@yahoo.com>
Date: Tue, 19 Feb 2013 08:23:11 +0000
Subject: [PATCH] PetscSolverNSCH added
---
AMDiS/src/parallel/PetscSolverNSCH.cc | 446 ++++++++++++++++++++++++++
AMDiS/src/parallel/PetscSolverNSCH.h | 243 ++++++++++++++
2 files changed, 689 insertions(+)
create mode 100644 AMDiS/src/parallel/PetscSolverNSCH.cc
create mode 100644 AMDiS/src/parallel/PetscSolverNSCH.h
diff --git a/AMDiS/src/parallel/PetscSolverNSCH.cc b/AMDiS/src/parallel/PetscSolverNSCH.cc
new file mode 100644
index 00000000..e97da33a
--- /dev/null
+++ b/AMDiS/src/parallel/PetscSolverNSCH.cc
@@ -0,0 +1,446 @@
+//
+// 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/PetscSolverNSCH.h"
+#include "parallel/PetscHelper.h"
+#include "TransformDOF.h"
+
+
+
+
+namespace AMDiS {
+
+ using namespace std;
+
+
+ PetscErrorCode pcShell(PC pc, Vec b, Vec x) // solve Px=b
+ {FUNCNAME("PCApply()");
+
+ void *ctx;
+ PCShellGetContext(pc, &ctx);
+ CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
+
+ /// extract vectors
+ Vec b1, b2, b34, b5, x1, x2, x34, x5, b3,b4;
+ data->globalMatrixSolver->extractVectorComponent(b, 0+3, &b1);
+ data->globalMatrixSolver->extractVectorComponent(b, 1+3, &b2);
+ data->globalMatrixSolver->extractVectorComponent(b, 0, &b34, 2);
+ data->globalMatrixSolver->extractVectorComponent(b, 2, &b5);
+
+ data->globalMatrixSolver->extractVectorComponent(x, 0+3, &x1);
+ data->globalMatrixSolver->extractVectorComponent(x, 1+3, &x2);
+ data->globalMatrixSolver->extractVectorComponent(x, 0, &x34, 2);
+ data->globalMatrixSolver->extractVectorComponent(x, 2, &x5);
+ data->globalMatrixSolver->extractVectorComponent(b, 0, &b3, 1);
+ data->globalMatrixSolver->extractVectorComponent(b, 1, &b4, 1);
+
+
+ /// solve Cahn-Hilliard Preconditioner
+ Vec y1, y2;
+ VecDuplicate(b1, &y1);
+ VecDuplicate(b2, &y2);
+
+ KSPSolve(data->kspMassCH, b1, y1); // M*y1 = b1
+ MatMultAdd(data->matMinusDeltaK, y1, b2, x1); // -> x1 := b2-delta*K*y1
+
+ KSPSolve(data->kspLaplaceCH, x1, y2); // (M+eps*sqrt(delta))*y2 = x1
+ MatMult(data->matMassCH, y2, x1); // x1 := M*y2
+
+ KSPSolve(data->kspLaplaceCH, x1, x2); // (M+eps*sqrt(delta))*x2 = x1
+ double factor = (*data->eps)/sqrt(*data->delta);
+ VecCopy(x2, x1); // x1 := x2
+ VecAXPBYPCZ(x1, 1.0, factor, -factor, y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1
+
+ VecDestroy(&y1);
+ VecDestroy(&y2);
+ /**/
+
+ /// solve Navier-Stokes Preconditioner
+ Vec tmp34, tmp5, tmp5_2, tmp34_2;
+ VecDuplicate(b34, &tmp34);
+ VecDuplicate(b34, &tmp34_2);
+ VecDuplicate(b5, &tmp5);
+ VecDuplicate(b5, &tmp5_2);
+
+ KSPSolve(data->kspVelocity, b34, tmp34);
+ VecScale(tmp34, -1.0);
+ MatMultAdd(data->matDiv, tmp34, b5, tmp5);
+
+ /// approximierte Schur Komplement
+ KSPSolve(data->kspLaplace, tmp5, x5);
+ { //project out constant Null-space
+ int vecSize;
+ VecGetSize(x5, &vecSize);
+ PetscScalar vecSum;
+ VecSum(x5, &vecSum);
+ vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
+ VecShift(x5, vecSum);
+ //VecView(y, PETSC_VIEWER_STDOUT_WORLD);
+ }
+ MatMult(data->matConDif, x5, tmp5);
+ KSPSolve(data->kspMass, tmp5, x5);
+ VecScale(x5,-1.0);
+
+ MatMultAdd(data->matGrad, x5, b34, tmp34);
+ KSPSolve(data->kspVelocity, tmp34, x34);
+ VecScale(x5,-1.0);
+/**/
+
+ VecDestroy(&tmp34);
+ VecDestroy(&tmp5);
+ VecDestroy(&tmp5_2);
+ VecDestroy(&tmp34_2);
+ VecDestroy(&b1);
+ VecDestroy(&b2);
+ VecDestroy(&b34);
+ VecDestroy(&b5);
+ VecDestroy(&x1);
+ VecDestroy(&x2);
+ VecDestroy(&x34);
+ VecDestroy(&x5);
+
+ }
+
+
+
+ PetscSolverNSCH::PetscSolverNSCH(string name)
+ : PetscSolverGlobalMatrix(name),
+ useOldInitialGuess(false),
+ laplaceSolutionMode(0),
+ massMatrixSolverCH(NULL),
+ laplaceMatrixSolverCH(NULL),
+ deltaKMatrixSolver(NULL),
+ pressureNullSpace(true),
+ velocitySolutionMode(0),
+ regularizeLaplace(0),
+ massSolutionMode(0),
+ massMatrixSolver(NULL),
+ laplaceMatrixSolver(NULL),
+ nu(NULL),
+ invTau(NULL),
+ solution(NULL),
+ phase(NULL)
+ {
+ Parameters::get(initFileStr + "->use old initial guess",
+ useOldInitialGuess);
+
+ Parameters::get(initFileStr + "->navierstokes->velocity solver",
+ velocitySolutionMode);
+
+ Parameters::get(initFileStr + "->navierstokes->mass solver",
+ massSolutionMode);
+
+ Parameters::get(initFileStr + "->navierstokes->laplace solver",
+ laplaceSolutionMode);
+ Parameters::get(initFileStr + "->navierstokes->regularize laplace",
+ regularizeLaplace);
+
+
+ }
+
+
+ void PetscSolverNSCH::solvePetscMatrix(SystemVector &vec,
+ AdaptInfo *adaptInfo)
+ {
+ FUNCNAME("PetscSolverNSCH::solvePetscMatrix()");
+
+ if (useOldInitialGuess) {
+ VecSet(getVecSolInterior(), 0.0);
+
+ for (int i = 0; i < solution->getSize(); i++)
+ setDofVector(getVecSolInterior(), solution->getDOFVector(i), i, true);
+
+ vecSolAssembly();
+ KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
+ }
+
+ PetscSolverGlobalMatrix::solvePetscMatrix(vec, adaptInfo);
+ }
+
+
+ void PetscSolverNSCH::initSolver(KSP &ksp)
+ {
+ FUNCNAME("PetscSolverNSCH::initSolver()");
+
+ // Create FGMRES based outer solver
+ MSG("CREATE POS 1: %p\n", &ksp);
+ KSPCreate(domainComm, &ksp);
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+ KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
+ petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCSHELL, 1e-6, 1e-8, 10000);
+ setConstantNullSpace(ksp, componentSpaces[0]->getMesh()->getDim() , true);
+ }
+
+
+ void PetscSolverNSCH::initPreconditioner(PC pc)
+ {
+ FUNCNAME("PetscSolverNSCH::initPreconditioner()");
+
+ MPI::COMM_WORLD.Barrier();
+ double wtime = MPI::Wtime();
+ int dim = componentSpaces[0]->getMesh()->getDim();
+ pressureComponent=dim;
+ const FiniteElemSpace *cahnHilliardFeSpace = componentSpaces[0+3];
+ const FiniteElemSpace *velocityFeSpace= componentSpaces[0];
+ const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
+
+ PCSetType(pc, PCSHELL);
+ PCShellSetApply(pc, pcShell);
+ PCShellSetContext(pc, &matShellContext);
+ matShellContext.globalMatrixSolver = (this);
+ matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
+
+
+ /// Init Cahn-Hilliard part
+
+ DOFMatrix laplaceMatrixCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
+ Operator laplaceOpCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
+
+ DOFMatrix massMatrixCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
+ Operator massOpCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
+
+ DOFMatrix deltaKMatrix(cahnHilliardFeSpace, cahnHilliardFeSpace);
+ Operator laplaceOp2(cahnHilliardFeSpace, cahnHilliardFeSpace);
+
+ {
+ Simple_ZOT zot;
+ massOpCH.addTerm(&zot);
+ laplaceOpCH.addTerm(&zot); // M
+ Simple_SOT sot2((*eps)*sqrt(*delta));
+ laplaceOpCH.addTerm(&sot2); // eps*sqrt(delta)*K
+ Simple_SOT sot(-(*delta));
+ laplaceOp2.addTerm(&sot); // -delta*K
+
+ massMatrixCH.assembleOperator(massOpCH);
+ massMatrixSolverCH = createSubSolver(0+3, "mass_");
+ massMatrixSolverCH->fillPetscMatrix(&massMatrixCH);
+ // === matrix (M + eps*sqrt(delta)*K) ===
+ laplaceMatrixCH.assembleOperator(laplaceOpCH);
+ laplaceMatrixSolverCH = createSubSolver(0+3, "laplace_");
+ laplaceMatrixSolverCH->fillPetscMatrix(&laplaceMatrixCH);
+
+ // === matrix (-delta*K) ===
+ deltaKMatrix.assembleOperator(laplaceOp2);
+ deltaKMatrixSolver = createSubSolver(0+3, "laplace2_");
+ deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
+ }
+
+ // === Setup solver ===
+ matShellContext.kspMassCH = massMatrixSolverCH->getSolver();
+ matShellContext.kspLaplaceCH = laplaceMatrixSolverCH->getSolver();
+ matShellContext.matMassCH = massMatrixSolverCH->getMatInterior();
+ matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();
+ matShellContext.eps = eps;
+ matShellContext.delta = delta;
+
+ petsc_helper::setSolver(matShellContext.kspMassCH, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
+ petsc_helper::setSolver(matShellContext.kspLaplaceCH, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+
+
+ /// Init Navier-Stokes part
+
+ DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
+ Operator massOp(pressureFeSpace, pressureFeSpace);
+ ZeroOrderTerm *massTerm = new Simple_ZOT;
+ massOp.addTerm(massTerm);
+ massMatrix.assembleOperator(massOp);
+ delete massTerm;
+ massMatrixSolver = createSubSolver(pressureComponent, "mass_");
+ massMatrixSolver->fillPetscMatrix(&massMatrix);
+ matShellContext.kspMass = massMatrixSolver->getSolver();
+
+ DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
+ Operator laplaceOp(pressureFeSpace, pressureFeSpace);
+ SecondOrderTerm *laplaceTerm = new Simple_SOT;
+ laplaceOp.addTerm(laplaceTerm);
+ 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));
+ vy.interpol(solution->getDOFVector(1));
+ if (dim == 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;
+ vp.interpol(solution->getDOFVector(dim+2));
+
+ densityFunctionTau = new LinearInterpolation(*rho1,*rho2,*invTau);
+ viscosityFunction = new LinearInterpolation(*nu1,*nu2);
+ densityFunction = new LinearInterpolation2(*rho1,*rho2);
+
+ conDif0 = new VecAtQP_ZOT(&vp, densityFunctionTau);
+ conDifOp.addTerm(conDif0);
+ conDif1 = new VecAtQP_SOT(&vp, viscosityFunction );
+ conDifOp.addTerm(conDif1);
+ conDif2 = new Vec2AtQP_FOT(&vx, &vp, densityFunction, 0);
+ conDifOp.addTerm(conDif2, GRD_PHI);
+
+ conDif3 = new Vec2AtQP_FOT(&vy, &vp, densityFunction, 1);
+ conDifOp.addTerm(conDif3, GRD_PHI);
+
+ if (dim == 3) {
+ conDif4 = new Vec2AtQP_FOT(&vz, &vp, densityFunction, 2);
+ conDifOp.addTerm(conDif4, GRD_PHI);
+ }
+
+ conDifMatrix.assembleOperator(conDifOp);
+
+ delete conDif0;
+ delete conDif1;
+ delete conDif2;
+ delete conDif3;
+ if (dim == 3) delete conDif4;
+ }
+
+ conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
+ conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
+ matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
+
+ extractMatrixComponent(mat[0][0], 2, 1, 0, 2, &(matShellContext.matDiv));
+ extractMatrixComponent(mat[0][0], 0, 2, 2, 1, &(matShellContext.matGrad));
+ extractMatrixComponent(mat[0][0], 0, 2, 0, 2, &(matShellContext.velocityMat));
+
+ ///erstelle kspVelocity
+ KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspVelocity));
+ KSPSetOperators(matShellContext.kspVelocity, matShellContext.velocityMat, matShellContext.velocityMat, SAME_NONZERO_PATTERN);
+
+ ///regularisiere LaplaceMatrix
+ if (regularizeLaplace)
+ {
+ PetscInt rows[1];
+ rows[0]=0;
+ MatZeroRows(laplaceMatrixSolver->getMatInterior(), 1, rows, 0, PETSC_NULL, PETSC_NULL);
+ KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspLaplace));
+ KSPSetOperators(matShellContext.kspLaplace, laplaceMatrixSolver->getMatInterior(), laplaceMatrixSolver->getMatInterior(), SAME_NONZERO_PATTERN);
+ }
+ else
+ { matShellContext.kspLaplace=laplaceMatrixSolver->getSolver();
+ setConstantNullSpace(matShellContext.kspLaplace);
+ }
+
+
+
+ // === Setup solver ===
+
+ 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, "mass_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
+ break;
+ default:
+ ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode);
+ }
+
+ switch (velocitySolutionMode) {
+ case 0:
+ petsc_helper::setSolver(matShellContext.kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+ break;
+ case 1:
+ petsc_helper::setSolverWithLu(matShellContext.kspVelocity, "", KSPPREONLY, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
+ break;
+ default:
+ ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
+ }
+
+
+ PCSetFromOptions(pc);
+
+ MSG("Setup of Cahn-Hilliard preconditioner needed %.5f seconds\n",
+ MPI::Wtime() - wtime);
+ }
+
+
+ void PetscSolverNSCH::exitPreconditioner(PC pc)
+ {
+ FUNCNAME("PetscSolverNSCH::exitPreconditioner()");
+
+ massMatrixSolverCH->destroyMatrixData();
+ massMatrixSolverCH->destroyVectorData();
+ laplaceMatrixSolverCH->destroyMatrixData();
+ laplaceMatrixSolverCH->destroyVectorData();
+ deltaKMatrixSolver->destroyMatrixData();
+ deltaKMatrixSolver->destroyVectorData();
+
+
+ massMatrixSolverCH->destroyVectorData();
+ laplaceMatrixSolverCH->destroyVectorData();
+ deltaKMatrixSolver->destroyVectorData();
+
+ delete massMatrixSolverCH;
+ massMatrixSolverCH = NULL;
+
+ delete laplaceMatrixSolverCH;
+ laplaceMatrixSolverCH = NULL;
+
+ delete deltaKMatrixSolver;
+ deltaKMatrixSolver = NULL;
+
+ 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;
+
+ KSPDestroy(&(matShellContext.kspVelocity));
+ if (regularizeLaplace)
+ KSPDestroy(&(matShellContext.kspLaplace));
+
+ MatDestroy(&(matShellContext.matGrad));
+ MatDestroy(&(matShellContext.matDiv));
+ MatDestroy(&(matShellContext.velocityMat));
+
+ delete densityFunction;
+ delete densityFunctionTau;
+ delete viscosityFunction;
+ }
+ }
+
diff --git a/AMDiS/src/parallel/PetscSolverNSCH.h b/AMDiS/src/parallel/PetscSolverNSCH.h
new file mode 100644
index 00000000..c5a0b926
--- /dev/null
+++ b/AMDiS/src/parallel/PetscSolverNSCH.h
@@ -0,0 +1,243 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// == http://www.amdis-fem.org ==
+// == ==
+// ============================================================================
+//
+// 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.
+
+
+
+/** \file PetscSolverNSCH.h */
+
+#include "parallel/PetscSolverGlobalMatrix.h"
+
+namespace AMDiS {
+
+ using namespace std;
+
+ struct CahnHilliardData2 {
+ KSP kspMassCH, kspLaplaceCH, kspVelocity, kspLaplace, kspMass;
+ Mat matMassCH, matMinusDeltaK, matGrad, matDiv, matConDif, matSchur, velocityMat;
+ PetscSolverGlobalMatrix *globalMatrixSolver;
+ double *eps, *delta;
+ MPI::Intracomm *mpiCommGlobal;
+ };
+
+
+ class PetscSolverNSCH : public PetscSolverGlobalMatrix
+ {
+ public:
+ /// Creator class
+
+ class IdFct : public AbstractFunction<double, double>
+ {
+ public:
+ IdFct()
+ : AbstractFunction<double, double>(0)
+ {}
+
+ double operator()(const double& x) const
+ {
+ return x;
+ }
+ };
+
+ class MultConstFct : public AbstractFunction<double, double>
+ {
+ public:
+ MultConstFct(double c)
+ : AbstractFunction<double, double>(1),
+ mConst(c)
+ {}
+
+ double operator()(const double& x) const
+ {
+ return mConst * x;
+ }
+ private:
+ double mConst;
+ };
+
+ class LinearInterpolation : public AbstractFunction<double, double>
+ {
+ public:
+ LinearInterpolation(double c1, double c2, double factor=1.0)
+ : AbstractFunction<double, double>(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 (result<cmin) result = cmin;
+ if (result>cmax) result = cmax;
+ return result;
+ }
+ private:
+ double a,b,cmin,cmax;
+ };
+
+
+ class LinearInterpolation2 : public BinaryAbstractFunction<double, double, double>
+ {
+ public:
+ LinearInterpolation2(double c1, double c2, double factor=1.0)
+ : BinaryAbstractFunction<double, double, double>(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 (result<cmin) result = cmin;
+ if (result>cmax) result = cmax;
+ return result * u;
+ }
+ private:
+ double a,b,cmin,cmax;
+ };
+
+
+ class Multiplier3 : public BinaryAbstractFunction<double, double, double>
+ {
+ public:
+ Multiplier3()
+ : BinaryAbstractFunction<double, double, double >(0)
+ {}
+
+ double operator()(const double& phi, const double& phase) const
+ {
+ return phase * phi;
+ }
+ };
+
+
+ class EinsMinus : public AbstractFunction<double, double>
+ {
+ public:
+ EinsMinus(double d)
+ : AbstractFunction<double, double>(0),
+ c(d)
+ {}
+
+ double operator()(const double& x) const
+ {
+ return c * std::max(1.0-x,0.000001);
+ }
+ private:
+ double c;
+ };
+
+ struct Multiplication : public BinaryAbstractFunction<double, double, double>{
+ double operator()(const double &v1, const double &v2) const{
+ return v2 / v1;
+ }
+ };
+
+
+
+
+ class Creator : public OEMSolverCreator
+ {
+ public:
+ virtual ~Creator() {}
+
+ /// Returns a new PetscSolver object.
+ OEMSolver* create()
+ {
+ return new PetscSolverNSCH(this->name);
+ }
+ };
+
+ PetscSolverNSCH(string name);
+
+ void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+
+ void setChData(double *epsPtr, double *deltaPtr)
+ {
+ eps = epsPtr;
+ delta = deltaPtr;
+ }
+
+
+ void setStokesData(double *invTauPtr, SystemVector *vec, double *nu1_=NULL, double *nu2_=NULL, double *rho1_=NULL, double *rho2_=NULL)
+ {
+ invTau = invTauPtr;
+ solution = vec;
+ nu1=nu1_;
+ nu2=nu2_;
+ rho1=rho1_;
+ rho2=rho2_;
+
+
+ }
+
+
+
+ void setPhase(DOFVector<double> *d, double eP3=0)
+ {
+ phase = d;
+ epsPhase3 = eP3;
+ }
+
+ protected:
+ void initSolver(KSP &ksp);
+
+ void initPreconditioner(PC pc);
+
+ void exitPreconditioner(PC pc);
+
+ private:
+ int pressureComponent, chComponent, velocityComponent;
+
+ 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;
+
+ /// 0: approximate solve 1: direct solver
+ int regularizeLaplace;
+
+ PetscSolver *massMatrixSolverCH, *laplaceMatrixSolverCH, *deltaKMatrixSolver;
+
+ CahnHilliardData2 matShellContext;
+
+ double *eps, *delta;
+
+ LinearInterpolation *densityFunctionTau;
+ LinearInterpolation *viscosityFunction;
+ LinearInterpolation2 *densityFunction;
+
+ PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver;
+
+ double *nu, *invTau, *nu1,*nu2,*rho1,*rho2;
+
+ double epsPhase3;
+
+ SystemVector* solution;
+
+ DOFVector<double>* phase;
+
+ IdFct idFct;
+
+ };
+
+}
+
--
GitLab