preconditioners added

extensions/preconditioner/CahnHilliard_.cc

+/******************************************************************************
+ *
+ * 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<double>);
+    }
+    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<double>);
+    }
+    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<double> 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<DOFAdmin*>(getMesh()->getVertexAdmin()),
+	Lagrange::getLagrange(getMesh()->getDim(), 1),
+	getMesh(),
+	"P1");
+      DOFVector<double> 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 >
+  
+}
+

extensions/preconditioner/CahnHilliard_.h

+/******************************************************************************
+ *
+ * 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<ProblemStat>
+{
+public: // definition of types
+
+  typedef BaseProblem<ProblemStat> 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

extensions/preconditioner/PetscSolverNavierStokes2.cc

+/******************************************************************************
+ *
+ * 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<NavierStokesSchurData2*>(ctx);
+
+    // Project out constant null space
+    {
+      int vecSize;
+      VecGetSize(x, &vecSize);
+      PetscScalar vecSum;
+      VecSum(x, &vecSum);
+      vecSum = vecSum / static_cast<PetscScalar>(-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<Matrix<DOFMatrix*> >& 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<double, WorldVector<double> >(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<int> velocityComponents;
+    for (size_t i = 0; i < static_cast<size_t>(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<double> vx(pressureFeSpace, "vx");
+    DOFVector<double> vy(pressureFeSpace, "vy");
+    DOFVector<double> vz(pressureFeSpace, "vz");
+    DOFVector<double> 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) {