diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 18155f135f3cd921f75cf825c7b2c873d66bc486..18245634d8332d18d3119ee1e9f24e3e8ac32fd1 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -509,7 +509,28 @@ namespace AMDiS {
PCFieldSplitSetIS(pc, splitName, is);
ISDestroy(&is);
}
+
+
+ void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::extractVectorComponent()");
+ IS is;
+ interiorMap->createIndexSet(is, i, numberOfComponents);
+ VecGetSubVector(input, is, output);
+ ISDestroy(&is);
+ }
+ void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::extractMatrixComponent()");
+ IS isrow, iscol;
+ interiorMap->createIndexSet(isrow, startRow, numberOfRows);
+ interiorMap->createIndexSet(iscol, startCol, numberOfCols);
+ MatGetSubMatrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
+ ISDestroy(&iscol);
+ ISDestroy(&isrow);
+ }
+
void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
{
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index 4fd425e339196a828dce4db6afbe06e0db101a47..dfde6305f772a18234b2a64728fd0e1976dd25f9 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -62,6 +62,10 @@ namespace AMDiS {
void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
void solveGlobal(Vec &rhs, Vec &sol);
+
+ void extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents=1);
+
+ void extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output);
void destroyMatrixData();
@@ -96,6 +100,8 @@ namespace AMDiS {
components.push_back(component);
createFieldSplit(pc, splitName, components);
}
+
+
virtual void initSolver(KSP &ksp);
diff --git a/extensions/base_problems/NavierStokesCahnHilliard.cc b/extensions/base_problems/NavierStokesCahnHilliard.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9001604b69d037396b93a1056608c930068e7e4e
--- /dev/null
+++ b/extensions/base_problems/NavierStokesCahnHilliard.cc
@@ -0,0 +1,418 @@
+#include "NavierStokesCahnHilliard.h"
+// #include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+
+using namespace AMDiS;
+
+NavierStokesCahnHilliard::NavierStokesCahnHilliard(const std::string &name_, bool createProblem) :
+ super(name_, createProblem),
+ useMobility(false),
+ doubleWell(0),
+ gamma(1.0),
+ eps(0.1),
+ minusEps(-0.1),
+ epsInv(10.0),
+ minusEpsInv(-10.0),
+ epsSqr(0.01),
+ minusEpsSqr(-0.01),
+ multiPhase(NULL),
+ densityPhase(NULL),
+ viscosityPhase(NULL),
+ viscosity1(1.0),
+ viscosity2(1.0),
+ density1(1.0),
+ density2(1.0),
+ refFunction(NULL),
+ refinement(NULL),
+ sigma(0.0),
+ surfaceTension(0.0)
+ {
+ dow = Global::getGeo(WORLD);
+ Initfile::get(name + "->viscosity1", viscosity1); // viscosity of fluid 1
+ Initfile::get(name + "->viscosity2", viscosity2); // viscosity of fluid 2
+ Initfile::get(name + "->density1", density1); // density of fluid 1
+ Initfile::get(name + "->density2", density2); // density of fluid 2
+
+ Initfile::get(name + "->non-linear term", nonLinTerm);
+ Initfile::get(name + "->theta", theta);
+ Initfile::get("adapt->timestep", tau);
+
+
+ Parameters::get(name + "->epsilon", eps); // interface width
+ Parameters::get(name + "->sigma", sigma);
+ double a_fac;
+ Parameters::get(name + "->a_factor", a_fac);
+ a=1.0/sigma*sqr(eps)*a_fac* (std::max(density1,density2)/tau);
+ b=1.0;
+ ab=a*b;
+ surfaceTension = -sigma/eps * 3.0/2.0/sqrt(2.0) * a ;
+
+ theta1 = 1.0-theta;
+ minusTheta1 = -theta1;
+
+ force.set(0.0);
+ Initfile::get(name + "->force", force);
+
+ // parameters for CH
+ Parameters::get(name + "->use mobility", useMobility); // mobility
+ Parameters::get(name + "->gamma", gamma); // mobility
+
+
+ // type of double well: 0= [0,1], 1= [-1,1]
+ Parameters::get(name + "->double-well type", doubleWell);
+
+ // transformation of the parameters
+ minusEps = -eps;
+ minus1 = -1.0;
+ epsInv = 1.0/eps;
+ minusEpsInv = -epsInv;
+ epsSqr = sqr(eps);
+ minusEpsSqr = -epsSqr*b;
+
+
+}
+
+
+NavierStokesCahnHilliard::~NavierStokesCahnHilliard()
+{ FUNCNAME("NavierStokesCahnHilliard::~NavierStokesCahnHilliard()");
+ delete densityPhase;
+ delete viscosityPhase;
+}
+
+
+void NavierStokesCahnHilliard::setTime(AdaptInfo* adaptInfo)
+ {
+ super::setTime(adaptInfo);
+ delta = gamma*adaptInfo->getTimestep();
+ }
+
+
+
+
+void NavierStokesCahnHilliard::initData()
+{ FUNCNAME("NavierStokes_TH_MultiPhase::initTimeInterface()");
+
+ #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ string initFileStr = name + "->space->solver", solverType = "";
+ Parameters::get(initFileStr, solverType);
+ if (solverType == "petsc-nsch") {
+ PetscSolverNSCH* solver = dynamic_cast<PetscSolverNSCH*>(prob->getSolver());
+ if (solver)
+ {
+ solver->setChData(&eps, &delta);
+ solver->setStokesData( getInvTau(), prob->getSolution(), &viscosity1, &viscosity2, &density1, &density2);
+ }
+ }
+ #endif
+
+ densityPhase = new DOFVector<double>(prob->getFeSpace(0), "densityPhase");
+ viscosityPhase = new DOFVector<double>(prob->getFeSpace(0), "viscosityPhase");
+
+ densityPhase->set(density1);
+ viscosityPhase->set(viscosity1);
+
+ if (velocity == NULL)
+ velocity = new DOFVector<WorldVector<double> >(getFeSpace(0), "velocity");
+
+ //fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace()->getMesh(), velocity);
+
+ super::initData();
+};
+
+
+void NavierStokesCahnHilliard::closeTimestep(AdaptInfo *adaptInfo)
+{ super::closeTimestep(adaptInfo);
+
+}
+
+void NavierStokesCahnHilliard::initTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("NavierStokes_TH_MultiPhase::initTimestep()");
+
+ super::initTimestep(adaptInfo);
+ refinement->refine(2);
+ LinearInterpolation1 dLI(density1, density2);
+ LinearInterpolation1 vLI(viscosity1, viscosity2);
+ transformDOF(prob->getSolution()->getDOFVector(dow+2), densityPhase, &dLI);
+ transformDOF(prob->getSolution()->getDOFVector(dow+2), viscosityPhase, &vLI);
+};
+
+
+
+void NavierStokesCahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo)
+ {
+ // meshFunction for klocal refinement around the interface of the phasefield-DOFVector
+ refFunction = new PhaseFieldRefinement(prob->getMesh());
+
+ if (getDoubleWellType() == 0) {
+ refinement = new RefinementLevelDOF(
+ prob->getFeSpace(),
+ refFunction,
+ prob->getSolution()->getDOFVector(dow+2)); // phaseField-DOFVector in [0,1]
+ } else {
+ refinement = new RefinementLevelDOF(
+ prob->getFeSpace(),
+ refFunction,
+ new PhaseDOFView<double>(prob->getSolution()->getDOFVector(dow+2))); // phaseField-DOFVector in [-1,1]
+ }
+
+ // initial refinement
+ refinement->refine(0);
+
+ for (int i = 0; i < 3; i++) {
+ solveInitialProblem2(adaptInfo); // update phaseField-DOFVector
+ refinement->refine((i < 4 ? 4 : 10)); // do k steps of local refinement
+ }
+
+ // solve all initial problems of the problems added to the CouplingTimeInterface
+ }
+
+
+
+void NavierStokesCahnHilliard::solveInitialProblem2(AdaptInfo *adaptInfo)
+{ FUNCNAME("NavierStokesCahnHilliard::solveInitialProblem()");
+
+ 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+3)->interpol(new Plane(a, dir));
+ }
+ else if (initialInterface == 1) {
+ /// schraege Linie
+ double theta = m_pi/4.0;
+ prob->getSolution()->getDOFVector(1+3)->interpol(new PlaneRotation(0.0, theta, 1.0));
+ transformDOFInterpolation(prob->getSolution()->getDOFVector(1+3),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+3)->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+3)->interpol(new Plane(a, -1.0));
+ transformDOFInterpolation(prob->getSolution()->getDOFVector(1+3),new Plane(b, 1.0), new AMDiS::Max<double>);
+ }
+ else if (initialInterface == 4) {
+ /// Kreis
+ double radius= 1.0, x=0, y=0;
+ Initfile::get(name + "->circle->radius", radius);
+ Initfile::get(name + "->circle->center_x", x);
+ Initfile::get(name + "->circle->center_y", y);
+ WorldVector<double> w;
+ w[0]=x; w[1]=y+adaptInfo->getTime();
+ prob->getSolution()->getDOFVector(1+3)->interpol(new Circle(radius, w));
+ } 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+3)->interpol(new Rectangle(width, height, center));
+ }
+
+ /// create phase-field from signed-dist-function
+ if (doubleWell == 0) {
+ transformDOF(prob->getSolution()->getDOFVector(1+3),
+ new SignedDistToPhaseField(initialEps));
+ } else {
+ transformDOF(prob->getSolution()->getDOFVector(1+3),
+ new SignedDistToCh(initialEps));
+ }
+ }
+}
+
+
+void NavierStokesCahnHilliard::fillOperators()
+{ FUNCNAME("NavierStokesCahnHilliard::fillOperators()");
+ MSG("NavierStokesCahnHilliard::fillOperators()\n");
+
+ const FiniteElemSpace* feSpace = prob->getFeSpace(0);
+
+ // c
+ Operator *opChMnew = new Operator(feSpace,feSpace);
+ opChMnew->addZeroOrderTerm(new Simple_ZOT);
+ Operator *opChMold = new Operator(feSpace,feSpace);
+ opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1+3)));
+ // -nabla*(grad(c))
+ Operator *opChL = new Operator(feSpace,feSpace);
+ opChL->addSecondOrderTerm(new Simple_SOT);
+
+ // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
+ Operator *opChLM = new Operator(feSpace,feSpace);
+ opChLM->addSecondOrderTerm(new Simple_SOT(gamma*ab));
+
+ // -2*c_old^3 + 3/2*c_old^2
+ Operator *opChMPowExpl = new Operator(feSpace,feSpace);
+ opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(1+3),
+ new Pow3Functor(-2.0)));
+ if (doubleWell == 0) {
+ opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(1+3),
+ new Pow2Functor(3.0/2.0)));
+ }
+
+ // -3*c_old^2 * c
+ Operator *opChMPowImpl = new Operator(feSpace,feSpace);
+ opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(1+3),
+ new Pow2Functor(-3.0)));
+ if (doubleWell == 0) {
+ opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getSolution()->getDOFVector(1+3),
+ new AMDiS::Factor<double>(3.0)));
+ opChMPowImpl->addZeroOrderTerm(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(*opChMnew,0+3,0+3, &ab); /// < phi*mu , psi >
+
+ prob->addMatrixOperator(*opChMPowImpl,0+3,1+3, &b); /// < -3*phi*c*c_old^2 , psi >
+ prob->addMatrixOperator(*opChL,0+3,1+3, &minusEpsSqr); /// < -eps^2*phi*grad(c) , grad(psi) >
+ // . . . vectorOperators . . . . . . . . . . . . . . .
+ prob->addVectorOperator(*opChMPowExpl,0+3, &b); /// < -2*phi*c_old^3 , psi >
+
+// setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
+
+ // dt(c) = laplace(mu) - u*grad(c)
+ // -----------------------------------
+ prob->addMatrixOperator(*opChMnew,1+3,1+3, &b); /// < phi*c/tau , psi >
+ prob->addMatrixOperator(*opChLM,1+3,0+3, getTau()); /// < phi*grad(mu) , grad(psi) >
+ // . . . vectorOperators . . . . . . . . . . . . . . .
+ prob->addVectorOperator(*opChMold,1+3, &b ); /// < phi*c^old/tau , psi >
+
+ /**/
+
+
+ /// Navier-Stokes part
+
+ WorldVector<DOFVector<double>* > vel;
+ for (unsigned i = 0; i < dow; i++){
+ vel[i]= prob->getSolution()->getDOFVector(i+2);
+ }
+
+
+ // fill operators for prob
+ for (unsigned i = 0; i < dow; ++i) {
+
+ /// < (1/tau)*rho*u'_i , psi >
+ Operator *opTime = new Operator(getFeSpace(i), getFeSpace(i));
+ if (density1==density2)
+ opTime->addTerm(new Simple_ZOT(density1));
+ else
+ opTime->addTerm(new VecAtQP_ZOT(densityPhase, NULL));
+ opTime->setUhOld(prob->getSolution()->getDOFVector(i));
+ prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());
+ prob->addVectorOperator(*opTime, i, getInvTau(), getInvTau());
+
+
+
+ /// < u^old*grad(u_i^old) , psi >
+/* Operator *opUGradU0 = new Operator(getFeSpace(i),getFeSpace(i));
+ if (density1==density2)
+ opUGradU0->addTerm(new WorldVec_FOT(vel, -1.0), GRD_PHI);
+ else
+ opUGradU0->addTerm(new WorldVecPhase_FOT(densityPhase, vel, -1.0), GRD_PHI);
+ opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
+ if (nonLinTerm == 0) {
+ prob->addVectorOperator(*opUGradU0, i);
+ } else {
+ prob->addVectorOperator(*opUGradU0, i, &theta1, &theta1);
+ }
+
+ if (nonLinTerm == 1) {
+ /// < u'*grad(u_i^old) , psi >
+ for (unsigned j = 2; j < 2+dow; ++j) {
+ Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(i));
+ if (density1==density2)
+ opUGradU1->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(i), j-2));
+ else
+ opUGradU1->addTerm(new VecAndPartialDerivative_ZOT( densityPhase, prob->getSolution()->getDOFVector(i), j-2));
+ prob->addMatrixOperator(*opUGradU1, i, j, &theta, &theta);
+ }
+ } else if (nonLinTerm == 2) {
+ /// < u^old*grad(u'_i) , psi >
+*/ for(unsigned j = 2; j < 2+dow; ++j) {
+ Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
+ opUGradU2->addTerm(new Vec2ProductPartial_FOT( densityPhase, prob->getSolution()->getDOFVector(j-2), j-2), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, i, i);
+ }
+
+ /**/
+ /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity
+ Operator *opLaplaceUi = new Operator(getFeSpace(i), getFeSpace(i));
+ opLaplaceUi->addTerm(new VecAtQP_SOT(viscosityPhase, NULL));
+ prob->addMatrixOperator(*opLaplaceUi, i, i);
+
+ /// < p , d_i(psi) >
+ Operator *opGradP = new Operator(getFeSpace(i),getFeSpace(2));
+ opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
+ prob->addMatrixOperator(*opGradP, i, 2);
+
+ /// external force, i.e. gravitational force
+ if (force[i]*force[i] > DBL_TOL) {
+ Operator *opForce = new Operator(getFeSpace(i), getFeSpace(i));
+ opForce->addZeroOrderTerm(new VecAtQP_ZOT(densityPhase, new Force(force[i])));
+ prob->addVectorOperator(*opForce, i);
+ }
+
+ /// < d_i(u'_i) , psi >
+ Operator *opDivU = new Operator(getFeSpace(2),getFeSpace(i));
+ opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+ prob->addMatrixOperator(*opDivU, 2, i);
+
+
+
+ ///coupling Operators
+ Operator *opNuGradC = new Operator(getFeSpace(i), getFeSpace(dow+1));
+ opNuGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(dow+2), i));
+ prob->addMatrixOperator(opNuGradC, i, dow+1, &surfaceTension, &surfaceTension);
+
+ Operator *opVGradC = new Operator(getFeSpace(dow+2), getFeSpace(i));
+ opVGradC->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(dow+2), i, b));
+ prob->addMatrixOperator(opVGradC, dow+2, i, getTau());
+ /**/
+
+ }
+
+ /**/
+};
+
+
+void NavierStokesCahnHilliard::addLaplaceTerm(int i)
+{ FUNCNAME("NavierStokes_TH_MultiPhase::addLaplaceTerm()");
+
+ /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
+ if (viscosity1!=viscosity2) {
+ for (unsigned j = 0; j < dow; ++j) {
+ Operator *opLaplaceUi1 = new Operator(getFeSpace(i), getFeSpace(j));
+ opLaplaceUi1->addTerm(new VecAtQP_IJ_SOT( viscosityPhase, NULL, j, i));
+ prob->addMatrixOperator(*opLaplaceUi1, i, j);
+ }
+ }/**/
+
+ /// < alpha*grad(u'_i) , grad(psi) >
+
+
+};
diff --git a/extensions/base_problems/NavierStokesCahnHilliard.h b/extensions/base_problems/NavierStokesCahnHilliard.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc879ec706e2dc08fb4143b170934cc8e8353b5f
--- /dev/null
+++ b/extensions/base_problems/NavierStokesCahnHilliard.h
@@ -0,0 +1,183 @@
+/** \file NavierStokesCahnHilliard.h */
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "POperators.h"
+#include "ExtendedProblemStat.h"
+#include "chns.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ #include "parallel/PetscSolverNSCH.h"
+#endif
+
+
+
+using namespace AMDiS;
+
+class NavierStokesCahnHilliard : public BaseProblem<ProblemStat>
+{
+public: // definition of types
+
+ typedef BaseProblem<ProblemStat> super;
+
+public: // public methods
+
+ NavierStokesCahnHilliard(const std::string &name_, bool createProblem = true);
+ ~NavierStokesCahnHilliard();
+
+ virtual void initData();
+
+ /** initTimestep of super and
+ * initialization of densityPhase and viscosityPhase
+ **/
+ virtual void initTimestep(AdaptInfo *adaptInfo);
+ virtual void closeTimestep(AdaptInfo *adaptInfo);
+
+ virtual void setTime(AdaptInfo* adaptInfo);
+
+/* * indicates the different phases. Will be initialized
+ * in \ref initTimeInterface with const 1
+ **/
+ void setMultiPhase(DOFVector<double>* p) { multiPhase=p; }
+
+
+ virtual void solveInitialProblem(AdaptInfo *adaptInfo);
+ virtual void solveInitialProblem2(AdaptInfo *adaptInfo);
+
+ double getEpsilon() { return eps; }
+ int getDoubleWellType() { return doubleWell; }
+
+protected: // protected methods
+
+ virtual void fillOperators();
+ virtual void addLaplaceTerm(int i);
+
+
+protected: // protected variables
+ ///viscosity of phase 1
+ double viscosity1;
+
+ ///viscosity of phase 2
+ double viscosity2;
+
+ ///density of phase 1
+ double density1;
+
+ ///density of phase 2
+ double density2;
+
+ double a, b, ab;
+
+ double tau;
+
+ /// phase dependent density
+ DOFVector<double> *densityPhase;
+
+ /// phase dependent viscosity
+ DOFVector<double> *viscosityPhase;
+
+ double delta;
+
+ /// phase inticator
+ DOFVector<double> *multiPhase;
+
+ DOFVector<WorldVector<double> >* velocity;
+
+ void calcVelocity()
+ {
+ if (dow == 1)
+ transformDOF(prob->getSolution()->getDOFVector(0), velocity, new AMDiS::Vec1WorldVec<double>);
+ else if (dow == 2)
+ transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), velocity, new AMDiS::Vec2WorldVec<double>);
+ else if (dow == 3)
+ transformDOF(prob->getSolution()->getDOFVector(0), prob->getSolution()->getDOFVector(1), prob->getSolution()->getDOFVector(2), velocity, new AMDiS::Vec3WorldVec<double>);
+ }
+
+ FileVectorWriter* fileWriter;
+
+ bool useMobility;
+
+ unsigned dow; // dimension of the world
+ unsigned dim;
+ PhaseFieldRefinement* refFunction;
+ RefinementLevelDOF *refinement;
+
+ double sigma, minus1; // coupling parameter to calculate the surface tension
+ double surfaceTension;// := sigma/epsilon
+
+ int doubleWell;
+
+ double gamma;
+ double eps;
+ double minusEps;
+ double epsInv;
+ double minusEpsInv;
+ double epsSqr;
+ double minusEpsSqr;
+
+ int nonLinTerm;
+ double theta;
+ double theta1;
+ double minusTheta1;
+ WorldVector<double> force;
+};
+
+
+
+class Force : public AbstractFunction<double, double>
+{
+public:
+ Force(double c_) :
+ AbstractFunction<double, double>(2), c(c_) {};
+
+ double operator()(const double &x) const {
+ return x*c;
+ }
+
+private:
+ double c;
+};
+
+
+class LinearInterpolation1 : public AbstractFunction<double, double>
+{
+public:
+
+ LinearInterpolation1(double c1, double c2) :
+ AbstractFunction<double, double>(1)
+ { a = (c1-c2)/2.0; b = (c1+c2)/2.0; cmin=std::min(c1,c2); cmax=std::max(c1,c2);}
+
+ double operator()(const double &phase) const {
+ double result = b+a*phase;
+ if (result<cmin) result = cmin;
+ if (result>cmax) result = cmax;
+ return result;
+ }
+
+private:
+ double a,b,cmin,cmax;
+};
+
+/** linear interpolation between two values (like density, viscosity)
+ * using a phase-field variable in [0,1]
+ **/
+class LinearInterpolation0 : public AbstractFunction<double, double>
+{
+public:
+
+ LinearInterpolation0(double val1_, double val2_) :
+ AbstractFunction<double, double>(1),
+ val1(val1_), val2(val2_) {}
+
+ double operator()(const double &phase) const {
+ return phase*val1 + (1.0-phase)*val2;
+ }
+
+private:
+
+ double val1;
+ double val2;
+};
+
+
diff --git a/extensions/demo/NavierStokesCahnHilliard_PC_coupled/CMakeLists.txt b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..438e779133fcbc7b93d2d7540ffa32946a0d6d80
--- /dev/null
+++ b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/CMakeLists.txt
@@ -0,0 +1,28 @@
+project("preconditioner")
+cmake_minimum_required(VERSION 2.6)
+unset(ENV{LIBRARY_PATH})
+find_package(AMDIS REQUIRED)
+if(AMDIS_FOUND)
+ message("AMDIS was found\n")
+ include(${AMDIS_USE_FILE})
+
+ SET(BASIS_LIBS ${AMDIS_LIBRARIES})
+endif(AMDIS_FOUND)
+
+# set(base_dir /home/spraetor/projects)
+#
+# file(GLOB common ${base_dir}/src/common/*.cc)
+# file(GLOB alglib ${base_dir}/src/extensions/alglib/cpp/src/*.cpp)
+#
+# include_directories(${base_dir}/src/common/)
+# include_directories(${meshconv_dir})
+
+
+ set(drivenCavity #src/NavierStokesCahnHilliard.cc
+ #src/PetscSolverNSCH.cc
+ src/benchmark.cc)
+ add_executable("benchmark" ${drivenCavity})
+ target_link_libraries("benchmark" ${BASIS_LIBS})
+
+
+
diff --git a/extensions/demo/NavierStokesCahnHilliard_PC_coupled/init/nsch.dat.2d b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/init/nsch.dat.2d
new file mode 100644
index 0000000000000000000000000000000000000000..07dabbd32a00082f0d3177c878896b45d19c84ed
--- /dev/null
+++ b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/init/nsch.dat.2d
@@ -0,0 +1,157 @@
+dimension of world: 2
+output_folder: .
+mesh_name: mesh
+
+% ====================== INCLUDES =========================
+#include "../init/reinit.inc.2d"
+
+
+% ============== USER-PARAMETER ==========================
+
+nsch->viscosity1: 1 % 1 0.1
+nsch->viscosity2: 10
+nsch->density1: 100 % 100 1
+nsch->density2: 1000
+nsch->sigma: 24.500 % 24.5 1.96
+nsch->a_factor: 1
+nsch->theta: 1
+nsch->force: [0.0, -0.98] % gravitational force [m/s^2]
+
+nsch->epsilon: 0.01
+nsch->initial epsilon: ${nsch->epsilon}
+nsch->gamma: 0.05
+nsch->transport term: 1.0
+nsch->initial interface: 4 % 0 für linie
+nsch->line->direction: 1
+nsch->line->pos: 1.0
+nsch->circle->radius: 0.25
+nsch->circle->center_x: 0.5
+nsch->circle->center_y: 0.5
+nsch->double-well type: 1
+nsch->use mobility: 0
+nsch->use conservation form: 0
+mesh->refinement->epsilon: ${nsch->epsilon}
+
+% ====================== TIMESTEPS ========================
+adapt->max iteration: 1
+adapt->max timestep iteration: 1
+adapt->max time iteration: 1
+
+adapt->timestep: 1.0e-3
+adapt->max timestep: 1e+10
+adapt->min timestep: 1e-6
+adapt->start time: 0.0
+adapt->end time: 3.0e0
+
+% ====================== MESH =============================
+
+mesh->refinement->initial level: 6 % 4 10
+mesh->refinement->level in inner domain: 6 % 4 10
+mesh->refinement->level in outer domain: 6 % 4 10
+mesh->refinement->level on interface: 13 % 8 14
+mesh->refinement->min inner interface value: 0.05
+mesh->refinement->max outer interface value: 0.95
+mesh->refinement->max inner interface value: 0.95
+mesh->refinement->min outer interface value: 0.05
+${mesh_name}->macro file name: ../macro/ns_ch.macro
+${mesh_name}->global refinements: 0
+${mesh_name}->check: 0
+nsch->space->mesh: ${mesh_name}
+
+% =========== OUTPUT ==============================================
+nsch->space->output[0]->filename: ${output_folder}/u1_
+nsch->space->output[1]->filename: ${output_folder}/u2_
+nsch->space->output[2]->filename: ${output_folder}/p_
+nsch->space->output[4]->filename: ${output_folder}/ch_
+
+
+% ============= PROBLEM-SPACES ==================================
+nsch->space->components: 5
+nsch->space->polynomial degree[0]: 2
+nsch->space->polynomial degree[1]: 2
+nsch->space->polynomial degree[2]: 1
+nsch->space->polynomial degree[3]: 2
+nsch->space->polynomial degree[4]: 2
+nsch->space->dim: 2
+
+% ================== SOLVER ======================================
+nsch->space->solver: petsc-nsch
+nsch->space->solver->navierstokes->regularize laplace: 1
+nsch->space->solver->use old initial guess: 1
+nsch->space->solver->navierstokes->velocity solver: 0
+nsch->space->solver->navierstokes->mass solver: 0
+nsch->space->solver->navierstokes->laplace solver: 0
+nsch->space->solver->use old initial guess: 1
+nsch->space->solver->symmetric strategy: 0
+nsch->space->solver->store symbolic: 0
+nsch->space->solver->ell: 8
+nsch->space->solver->max iteration: 200
+nsch->space->solver->tolerance: 1.e-10
+nsch->space->solver->info: 1
+nsch->space->solver->left precon: ilu
+
+%nsch->space->solver: umfpack
+%nsch->space->solver->symmetric strategy: 0
+%nsch->space->solver->store symbolic: 0
+%nsch->space->solver->ell: 8
+%nsch->space->solver->max iteration: 200
+%nsch->space->solver->tolerance: 1.e-8
+%nsch->space->solver->info: 1
+%nsch->space->solver->left precon: ilu
+
+% =================== OUTPUT =========================================
+
+nsch->space->output[0]->ParaView animation: 1
+nsch->space->output[0]->ParaView format: 1
+nsch->space->output[0]->write every i-th timestep: 1
+%nsch->space->output->compression: gzip
+nsch->space->output[0]->append index: 1
+nsch->space->output[0]->index length: 9
+nsch->space->output[0]->index decimals: 7
+
+nsch->space->output[1]->ParaView animation: 1
+nsch->space->output[1]->ParaView format: 1
+nsch->space->output[1]->write every i-th timestep: 1
+%nsch->space->output->compression: gzip
+nsch->space->output[1]->append index: 1
+nsch->space->output[1]->index length: 9
+nsch->space->output[1]->index decimals: 7
+
+nsch->space->output[2]->ParaView animation: 1
+nsch->space->output[2]->ParaView format: 1
+nsch->space->output[2]->write every i-th timestep: 1
+%nsch->space->output->compression: gzip
+nsch->space->output[2]->append index: 1
+nsch->space->output[2]->index length: 9
+nsch->space->output[2]->index decimals: 7
+
+nsch->space->output[4]->ParaView animation: 1
+nsch->space->output[4]->ParaView format: 1
+nsch->space->output[4]->write every i-th timestep: 1
+%nsch->space->output->compression: gzip
+nsch->space->output[4]->append index: 1
+nsch->space->output[4]->index length: 9
+nsch->space->output[4]->index decimals: 7
+
+
+
+
+nsch->space->output->write serialization: 0
+nsch->space->output->serialization filename: serial_ch
+nsch->space->input->read serialization: 0
+nsch->space->input->serialization filename: serial_ch
+
+
+% ====================== MAIN INITFILE ====================
+
+% helper problem to initialize all FeSpaces
+nsch->mesh: ${mesh_name}
+nsch->dim: 2
+
+
+
+% ====================== ESTIMATORS =======================
+adapt->strategy: 0 % 0=explicit, 1=implicit
+
+WAIT: 1
+
diff --git a/extensions/demo/NavierStokesCahnHilliard_PC_coupled/init/reinit.inc.2d b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/init/reinit.inc.2d
new file mode 100644
index 0000000000000000000000000000000000000000..4d7d822ea427b31fcc1ae530635259e800a8e4ca
--- /dev/null
+++ b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/init/reinit.inc.2d
@@ -0,0 +1,5 @@
+reinit->tolerance: 1.e-4
+reinit->maximal number of iteration steps: 100
+reinit->Gauss-Seidel iteration: 1
+reinit->infinity value: 1.e8
+reinit->boundary initialization: 3
diff --git a/extensions/demo/NavierStokesCahnHilliard_PC_coupled/macro/ns_ch.macro b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/macro/ns_ch.macro
new file mode 100644
index 0000000000000000000000000000000000000000..8fe05a720ae77bf1a2f22a90c2d20e3db44edf95
--- /dev/null
+++ b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/macro/ns_ch.macro
@@ -0,0 +1,37 @@
+DIM: 2
+DIM_OF_WORLD: 2
+
+number of elements: 8
+number of vertices: 8
+
+element vertices:
+0 1 6
+1 2 7
+2 3 7
+3 4 7
+4 5 6
+5 0 6
+1 4 6
+4 1 7
+
+
+element boundaries:
+0 0 1
+0 0 1
+0 0 2
+0 0 1
+0 0 1
+0 0 2
+0 0 0
+0 0 0
+
+
+vertex coordinates:
+ 1.0 0.0
+ 1.0 1.0
+ 1.0 2.0
+ 0.0 2.0
+ 0.0 1.0
+ 0.0 0.0
+ 0.5 0.5
+ 0.5 1.5
diff --git a/extensions/demo/NavierStokesCahnHilliard_PC_coupled/run b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/run
new file mode 100644
index 0000000000000000000000000000000000000000..2fbc804b20e5a487fd0c02e7e36c544d38cef674
--- /dev/null
+++ b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/run
@@ -0,0 +1,19 @@
+# This demo implements the Navier-Stokes-Cahn-Hilliard Benchmark from
+# S. Aland, A. Voigt. Benchmark computations of diffuse-interface models for two-dimensional bubble dynamics. Int. J. Num. Meth. Fluids (2012)
+
+output="output_parallel"
+initfile="init/nsch.dat.2d"
+mkdir $output
+cd $output
+mkdir serials
+cp -r ../src .
+cp ../$initfile .
+
+mpiexec -n 2 ../benchmark ../$initfile -ns_ksp_atol 1e-7 -ns_ksp_rtol 0 -ch_ksp_atol 1e-8 -ch_ksp_rtol 0 -laplace_pc_type hypre -laplace_pc_hypre_boomeramg_relax_type_coarse symmetric-SOR/Jacobi
+
+# sequentiell (solver noch auf umfpack stellen) -> #../nsch ../$initfile
+
+
+###### zum parallel debugging
+###### mpiexec -n 2 valgrind --tool=memcheck -q --num-callers=20 --log-file=valgrind.log.%p ../drivenCavity ../$initfile -malloc off -log_summary
+
diff --git a/extensions/demo/NavierStokesCahnHilliard_PC_coupled/src/benchmark.cc b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/src/benchmark.cc
new file mode 100644
index 0000000000000000000000000000000000000000..503eebb6345753561059a1598e14ec42d6051707
--- /dev/null
+++ b/extensions/demo/NavierStokesCahnHilliard_PC_coupled/src/benchmark.cc
@@ -0,0 +1,69 @@
+#include "AMDiS.h"
+#include "NavierStokesCahnHilliard.h"
+//#include "NavierStokes_TH_MultiPhase.h"
+//#include "CahnHilliardNavierStokes_.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+
+
+class MyNavierStokesCahnHilliard : public NavierStokesCahnHilliard
+{
+public:
+ MyNavierStokesCahnHilliard(std::string name) : NavierStokesCahnHilliard(name)
+ { }
+
+ void fillBoundaryConditions()
+ { FUNCNAME("MyNavierStokes::fillBoundaryConditions()");
+
+ DOFVector<double> *zeroDOF = new DOFVector<double>(getFeSpace(0), "zero");
+ zeroDOF->set(0.0);
+ size_t dow = Global::getGeo(WORLD);
+
+ /// at rigid wall: no-slip boundary condition
+ // oben-unten
+ getProblem()->addDirichletBC(2, 0, 0, zeroDOF);
+ getProblem()->addDirichletBC(2, 1, 1, zeroDOF);
+ // links-rechts
+ getProblem()->addDirichletBC(1, 0, 0, zeroDOF);
+ //getProblem()->addDirichletBC(1, 1, 1, zeroDOF);
+ }
+
+};
+
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+ AMDiS::init(argc, argv);
+
+ #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ CreatorMap<OEMSolver>::addCreator("petsc-nsch", new PetscSolverNSCH::Creator);
+ #endif
+
+ MyNavierStokesCahnHilliard prob("nsch");
+
+ // Adapt-Infos
+ AdaptInfo adaptInfo("adapt", prob.getNumComponents());
+ prob.initialize(INIT_ALL);
+ prob.initData();
+ prob.initTimeInterface();
+
+ AdaptInstationary adaptInstat("adapt", prob, adaptInfo, prob, adaptInfo);
+
+ ptime start_time = microsec_clock::local_time();
+ int error_code = adaptInstat.adapt();
+ time_duration td = microsec_clock::local_time()-start_time;
+
+ MSG("elapsed time= %d sec\n", td.total_seconds());
+
+ AMDiS::finalize();
+
+ return error_code;
+};