blub

AMDiS/src/io/VtkWriter.cc

@@ -161,6 +161,8 @@ namespace AMDiS {
   {
     FUNCNAME("VtkWriter::writeFile()");
 
+    return;
+
     DataCollector<> dc(values->getFeSpace(), values);
     vector<DataCollector<>*> dcList(0);
     dcList.push_back(&dc);

AMDiS/src/parallel/PetscHelper.cc

@@ -266,6 +266,7 @@ namespace AMDiS {
       PC pc;
       KSPGetPC(ksp, &pc);
       PCSetType(pc, pcType);
+      PCSetFromOptions(pc);
     }
   }
 

AMDiS/src/parallel/PetscProblemStat.cc

@@ -50,7 +50,7 @@ namespace AMDiS {
 #else
       ERROR_EXIT("AMDiS was compiled without BDDC-ML support!\n");
 #endif
-    } else if (tmp == "petsc-stokes") {
+    } else if (tmp == "petsc-navierstokes") {
       petscSolver = new PetscSolverNavierStokes(initFileStr);
     } else {
       ERROR_EXIT("No parallel solver %s available!\n", tmp.c_str());

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -841,6 +841,8 @@ namespace AMDiS {
 						     int constFeSpace,
 						     bool test)
   {
+    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");
+
     Vec nullSpaceBasis;
     VecDuplicate(getVecSolInterior(), &nullSpaceBasis);
     
@@ -871,4 +873,14 @@ namespace AMDiS {
   }
 
 
+  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
+  {
+    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");
+
+    MatNullSpace matNullSpace;
+    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
+    KSPSetNullSpace(ksp, matNullSpace);
+    MatNullSpaceDestroy(&matNullSpace);
+  }
+
 }

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -123,6 +123,8 @@ namespace AMDiS {
 
     void setConstantNullSpace(KSP ksp, int constFeSpace, bool test = false);
 
+    void setConstantNullSpace(KSP ksp);
+
   protected:
     bool zeroStartVector;
 

AMDiS/src/parallel/PetscSolverNavierStokes.cc

@@ -13,6 +13,7 @@
 
 #include "parallel/PetscSolverNavierStokes.h"
 #include "parallel/PetscHelper.h"
+#include "TransformDOF.h"
 
 namespace AMDiS {
 
@@ -38,9 +39,10 @@ namespace AMDiS {
       laplaceMatrixSolver(NULL),
       nu(NULL),
       invTau(NULL),
-      solution(NULL)
+      solution(NULL),
+      phase(NULL)
   {
-    Parameters::get(initFileStr + "->stokes->pressure component", 
+    Parameters::get(initFileStr + "->navierstokes->pressure component", 
 		    pressureComponent);
 
     TEST_EXIT(pressureComponent >= 0)
@@ -61,7 +63,7 @@ namespace AMDiS {
     petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
     
     // Create null space information.
-    setConstantNullSpace(ksp, pressureComponent);
+    setConstantNullSpace(ksp, pressureComponent, true);
   }
 
 
@@ -69,6 +71,10 @@ namespace AMDiS {
   {
     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<int> velocityComponents;
     velocityComponents.push_back(0);
     velocityComponents.push_back(1);
@@ -92,7 +98,7 @@ namespace AMDiS {
     KSP kspSchur = subKsp[1];
     PetscFree(subKsp);
 
-    petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+    petsc_helper::setSolver(kspVelocity, "velocity_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
 
     KSPSetType(kspSchur, KSPPREONLY);
     PC pcSub;
@@ -101,10 +107,7 @@ namespace AMDiS {
     PCShellSetApply(pcSub, pcSchurShell);
     PCShellSetContext(pcSub, &matShellContext);
 
-    MatNullSpace matNullSpace;
-    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
-    KSPSetNullSpace(kspSchur, matNullSpace);
-    MatNullSpaceDestroy(&matNullSpace);
+    setConstantNullSpace(kspSchur);
 
 
     // === Mass matrix solver ===
@@ -112,19 +115,25 @@ namespace AMDiS {
     const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
     DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
     Operator massOp(pressureFeSpace, pressureFeSpace);
-    Simple_ZOT zot;
-    massOp.addTerm(&zot);
+    //    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);
-    Simple_SOT sot;
-    laplaceOp.addTerm(&sot);
+    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);
@@ -134,19 +143,39 @@ namespace AMDiS {
 
     DOFVector<double> vx(pressureFeSpace, "vx");
     DOFVector<double> vy(pressureFeSpace, "vy");
+    DOFVector<double> vp(pressureFeSpace, "vp");
     vx.interpol(solution->getDOFVector(0));
     vy.interpol(solution->getDOFVector(1));
 
     DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
     Operator conDifOp(pressureFeSpace, pressureFeSpace);
-    Simple_ZOT conDif0(*invTau);
-    conDifOp.addTerm(&conDif0);
-    Simple_SOT conDif1(*nu);
-    conDifOp.addTerm(&conDif1);
-    VecAtQP_FOT conDif2(&vx, &idFct, 0);
-    conDifOp.addTerm(&conDif2, GRD_PHI);
-    VecAtQP_FOT conDif3(&vy, &idFct, 1);
-    conDifOp.addTerm(&conDif3, GRD_PHI);
+
+    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);
 
@@ -160,8 +189,30 @@ namespace AMDiS {
     matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
     matShellContext.matConDif = conDifMatrixSolver->getMatInterior();    
 
-    petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
-    petsc_helper::setSolver(matShellContext.kspLaplace, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+    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);
+
   }
 
 

AMDiS/src/parallel/PetscSolverNavierStokes.h

@@ -42,7 +42,7 @@ namespace AMDiS {
     {
     public:
       IdFct() 
-	: AbstractFunction<double, double>(0)
+	: AbstractFunction<double, double>(1)
       {}
       
       double operator()(const double& x) const 
@@ -51,6 +51,60 @@ namespace AMDiS {
       }
     };
 
+    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 Multiplier3 : public BinaryAbstractFunction<double, double, double>
+    {
+    public:
+      Multiplier3() 
+	: BinaryAbstractFunction<double, double, double >(6) 
+      {}
+      
+      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>(2),
+	  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;
+      }
+    };
+    
+
+
   public:
     PetscSolverNavierStokes(string name);
 
@@ -61,6 +115,11 @@ namespace AMDiS {
       solution = vec;
     }
     
+    void setPhase(DOFVector<double> *d)
+    {
+      phase = d;
+    }
+
   protected:
     void initSolver(KSP &ksp);
 
@@ -77,6 +136,8 @@ namespace AMDiS {
 
     SystemVector* solution;
 
+    DOFVector<double>* phase;
+    
     IdFct idFct;
   };
 

AMDiS/src/reinit/HL_SignedDistTraverse.cc

@@ -80,12 +80,12 @@ void HL_SignedDistTraverse::initializeBoundary()
 	for (int i = 0; i <= dim; i++) {
 	  
 	  // --> for test purposes:
-	  if (distVec[i] > 1000)
-	    cout << "\nElement: " << elInfo->getElement()->getIndex() 
-		 << ", Knoten: " << i << " , keine Randwertinitialisierung !\n";
+	  if (distVec[i] > 1000) {
+	    MSG("Element %d  Knoten %d  keine Randwertinitialisierung!\n",
+		elInfo->getElement()->getIndex(), i);
+	  }
 		  
 	  // --> end: for test purposes	  
-	  
 	  if ((*sD_DOF)[locInd[i]] > distVec[i]) {
 	    (*sD_DOF)[locInd[i]] = distVec[i];
 	    //If Distance is corrected, calculate new velocity.
@@ -103,6 +103,8 @@ void HL_SignedDistTraverse::initializeBoundary()
 
 void HL_SignedDistTraverse::HL_updateIteration()
 {
+  FUNCNAME("HL_SignedDistTraverse::HL_updateIteration()");
+
   // ===== Create DOF vector for the last iteration step. =====
   if (sDOld_DOF)
     delete sDOld_DOF;
@@ -140,14 +142,14 @@ void HL_SignedDistTraverse::HL_updateIteration()
     sDOld_DOF->copy(const_cast<DOFVector<double> &>(*sD_DOF));
   }
   
-  cout << "\n\nCalculation of signed distance function via mesh traverse iteration:\n";
 
+  MSG("Calculation of signed distance function via mesh traverse iteration:\n");
   if (GaussSeidelFlag)
-    cout << "\tGauss-Seidel iteration\n";  
+    MSG("\tGauss-Seidel iteration\n");
   else
-    cout << "\tJacobi iteration\n";
+    MSG("\tJacobi iteration\n");
   
-  cout << "\tnumber of iterations needed: " << itCntr << "\n\n";
+  MSG("\tnumber of iterations needed: %d\n", itCntr);
 }