Commit ff13277f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

correction of MeshManipulation

parent d2ebd11b
...@@ -46,9 +46,9 @@ public: ...@@ -46,9 +46,9 @@ public:
~CouplingBaseProblem() { } ~CouplingBaseProblem() { }
virtual void initialize(Flag initFlag, void initialize(Flag initFlag,
ProblemStat *adoptProblem = NULL, ProblemStatSeq *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING) Flag adoptFlag = INIT_NOTHING) _OVERRIDE_
{ {
for (size_t i = 0; i < baseProblems.size(); i++) { for (size_t i = 0; i < baseProblems.size(); i++) {
for (size_t j = 0; j < baseProblems[i]->getNumProblems(); j++) for (size_t j = 0; j < baseProblems[i]->getNumProblems(); j++)
......
...@@ -83,6 +83,8 @@ public: ...@@ -83,6 +83,8 @@ public:
#define BASEPROBLEM_ARGS(z, n, text) BOOST_PP_COMMA_IF(n) _BaseProblemType ## n &prob ## n #define BASEPROBLEM_ARGS(z, n, text) BOOST_PP_COMMA_IF(n) _BaseProblemType ## n &prob ## n
#define BASEPROBLEM_TYPES(z, n, text) class _BaseProblemType ## n #define BASEPROBLEM_TYPES(z, n, text) class _BaseProblemType ## n
// constructors created by C-Macros for each number of arguments, using BOOST_PP-Macros
#define COUPLED_BASEPROBLEM_CONSTRUCTOR(z, n, text) \ #define COUPLED_BASEPROBLEM_CONSTRUCTOR(z, n, text) \
template<BOOST_PP_ENUM(BOOST_PP_INC(n), BASEPROBLEM_TYPES, BaseProblemType)> \ template<BOOST_PP_ENUM(BOOST_PP_INC(n), BASEPROBLEM_TYPES, BaseProblemType)> \
...@@ -110,7 +112,7 @@ public: ...@@ -110,7 +112,7 @@ public:
~CouplingBaseProblem() { } ~CouplingBaseProblem() { }
void initialize(Flag initFlag, void initialize(Flag initFlag,
ProblemStat *adoptProblem = NULL, ProblemStatSeq *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING) _OVERRIDE_ Flag adoptFlag = INIT_NOTHING) _OVERRIDE_
{ {
tools::FOR_EACH< detail::AddProblem >::loop2(*this, baseProblems); tools::FOR_EACH< detail::AddProblem >::loop2(*this, baseProblems);
......
...@@ -88,9 +88,18 @@ public: ...@@ -88,9 +88,18 @@ public:
} }
~CouplingBaseProblem() { } ~CouplingBaseProblem() { }
/**
* Add the problems to the iterationInterface, timeInterface and couplingProblemStat.
* As a consequence all problem can be initialized one after another and in the
* adaption loop they are solved in rotation.
*
* In the adaption loop the problems are solved the same order as they are added to the
* iterationInterface in this method. This order can be changed manually in the oneIteration
* method.
**/
void initialize(Flag initFlag, void initialize(Flag initFlag,
ProblemStat *adoptProblem = NULL, ProblemStatSeq *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING) _OVERRIDE_ Flag adoptFlag = INIT_NOTHING) _OVERRIDE_
{ {
tools::FOR_EACH< detail::AddProblem >::loop2(*this, baseProblems); tools::FOR_EACH< detail::AddProblem >::loop2(*this, baseProblems);
...@@ -106,6 +115,12 @@ public: ...@@ -106,6 +115,12 @@ public:
virtual void initData() {} virtual void initData() {}
/**
* At first the initData method is called for all baseProblems, then
* the problems are filled with operators and coupling operators as well as
* boundary conditions are added.
**/
virtual void initTimeInterface() virtual void initTimeInterface()
{ {
tools::FOR_EACH< detail::InitData >::loop(baseProblems); tools::FOR_EACH< detail::InitData >::loop(baseProblems);
......
...@@ -28,7 +28,7 @@ public: ...@@ -28,7 +28,7 @@ public:
typedef BaseProblem_RB BaseProblemType; typedef BaseProblem_RB BaseProblemType;
typedef CouplingProblemStatImpl<ExtendedRosenbrockStationary> CProblemStat; typedef CouplingProblemStatImpl<ExtendedRosenbrockStationary> CProblemStat;
CoupledBaseProblem_RB(std::string name_, CouplingBaseProblem_RB(std::string name_,
BaseProblemType *prob0_, BaseProblemType *prob0_,
BaseProblemType *prob1_=NULL, BaseProblemType *prob1_=NULL,
BaseProblemType *prob2_=NULL, BaseProblemType *prob2_=NULL,
...@@ -43,11 +43,11 @@ public: ...@@ -43,11 +43,11 @@ public:
} }
} }
~CoupledBaseProblem_RB() { } ~CouplingBaseProblem_RB() { }
/// get the j-th stage-solution-vector of the i-th problem /// get the j-th stage-solution-vector of the i-th problem
DOFVector<double> *getStageSolution(int i, int j) DOFVector<double> *getStageSolution(int i, int j)
{ FUNCNAME("CoupledBaseProblem::getStageSolution(i,j)"); { FUNCNAME("CouplingBaseProblem_RB::getStageSolution(i,j)");
TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents()) TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents())
("Indices out of range!\n"); ("Indices out of range!\n");
...@@ -56,7 +56,7 @@ public: ...@@ -56,7 +56,7 @@ public:
/// get the j-th unVec-vector of the i-th problem /// get the j-th unVec-vector of the i-th problem
DOFVector<double> *getUnVec(int i, int j=0) DOFVector<double> *getUnVec(int i, int j=0)
{ FUNCNAME("CoupledBaseProblem::getUnVec(i,j)"); { FUNCNAME("CouplingBaseProblem_RB::getUnVec(i,j)");
TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents()) TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents())
("Indices out of range!\n"); ("Indices out of range!\n");
......
...@@ -31,6 +31,7 @@ ch->space->dim: 2 ...@@ -31,6 +31,7 @@ ch->space->dim: 2
% ================== SOLVER ====================================== % ================== SOLVER ======================================
ch->space->solver: petsc-ch % umfpack %petsc-ch ch->space->solver: petsc-ch % umfpack %petsc-ch
ch->space->solver->backend: no
ch->space->solver->use old initial guess: 1 ch->space->solver->use old initial guess: 1
ch->space->solver->symmetric strategy: 0 ch->space->solver->symmetric strategy: 0
ch->space->solver->store symbolic: 0 ch->space->solver->store symbolic: 0
......
...@@ -48,6 +48,7 @@ ns->space->dim: 2 ...@@ -48,6 +48,7 @@ ns->space->dim: 2
% ================== SOLVER ====================================== % ================== SOLVER ======================================
ns->space->solver: petsc-navierstokes %umfpack % petsc-navierstokes ns->space->solver: petsc-navierstokes %umfpack % petsc-navierstokes
ns->space->solver->backend: no
ns->space->solver->navierstokes->use old initial guess: 1 ns->space->solver->navierstokes->use old initial guess: 1
parallel->use zero start vector: 0 parallel->use zero start vector: 0
ns->space->solver->navierstokes->pressure component: 2 ns->space->solver->navierstokes->pressure component: 2
......
...@@ -18,10 +18,14 @@ ...@@ -18,10 +18,14 @@
#ifndef CAHN_HILLIARD_NAVIER_STOKES_H #ifndef CAHN_HILLIARD_NAVIER_STOKES_H
#define CAHN_HILLIARD_NAVIER_STOKES_H #define CAHN_HILLIARD_NAVIER_STOKES_H
#include "AMDiS.h"
// coupling structures // coupling structures
#include "CouplingIterationInterface.h" #if __cplusplus > 199711L
#include "CouplingTimeInterface.h" #include "CouplingBaseProblem2_cxx11.h"
#include "CouplingProblemStat.h" #else
#include "CouplingBaseProblem2.h"
#endif
#include "CahnHilliard_.h" #include "CahnHilliard_.h"
...@@ -33,38 +37,6 @@ ...@@ -33,38 +37,6 @@
using namespace AMDiS; using namespace AMDiS;
class ID : public AbstractFunction<double, double>
{
public:
ID()
: AbstractFunction<double, double>(2) {};
double operator()(const double& x) const {
return x;
}
};
class Y : public AbstractFunction<double, WorldVector<double> >
{
public:
Y()
: AbstractFunction<double, WorldVector<double> >(1) {};
double operator()(const WorldVector<double>& x) const {
return x[1];
}
};
/** /**
* \ingroup Problem * \ingroup Problem
* *
...@@ -72,90 +44,87 @@ public: ...@@ -72,90 +44,87 @@ public:
*/ */
template<typename CH_Type = CahnHilliard_, template<typename CH_Type = CahnHilliard_,
typename NS_Type = NavierStokes_TaylorHood> typename NS_Type = NavierStokes_TaylorHood>
class CahnHilliardNavierStokes_ : public CouplingIterationInterface, class CahnHilliardNavierStokes_ : public CouplingBaseProblem<ProblemStat, CH_Type, NS_Type>
public CouplingTimeInterface,
public CouplingProblemStat
{ {
public:
typedef CouplingBaseProblem<ProblemStat, CH_Type, NS_Type> super;
public: public:
CahnHilliardNavierStokes_(std::string name_, CH_Type *chProb_, NS_Type *nsProb_) CahnHilliardNavierStokes_(std::string name_, CH_Type *chProb_, NS_Type *nsProb_)
: CouplingProblemStat(name_), : super(name_, *chProb_, *nsProb_),
chProb(chProb_), chProb(chProb_),
nsProb(nsProb_), nsProb(nsProb_),
refFunction(NULL), refFunction(NULL),
refinement(NULL), refinement(NULL),
sigma(0.0), sigma(0.0),
surfaceTension(0.0) surfaceTension(0.0),
y(NULL),
yc1(NULL),
yc2(NULL)
{ {
dow = Global::getGeo(WORLD); Parameters::get(super::name + "->sigma", sigma);
Parameters::get(name + "->sigma", sigma);
surfaceTension = sigma/chProb->getEpsilon() * 3.0/2.0/sqrt(2.0); surfaceTension = sigma/chProb->getEpsilon() * 3.0/2.0/sqrt(2.0);
} }
~CahnHilliardNavierStokes_() { ~CahnHilliardNavierStokes_()
if (refFunction != NULL) {
if (refFunction) {
delete refFunction; delete refFunction;
if (refinement != NULL) refFunction = NULL;
}
if (refinement) {
delete refinement; delete refinement;
refinement = NULL;
}
if (y) {
delete y;
y = NULL;
}
if (yc1) {
delete yc1;
yc1 = NULL;
}
if (yc2) {
delete yc2;
yc2 = NULL;
}
} }
/** /**
* Add the problems to the iterationInterface, timeInterface and couplingProblemStat. * All needed DOFVectors are initialized so that they can be used inside
* As a consequence all problem can be initialized one after another and in the * of the operators.
* adaption loop they are solved in rotation.
* At the end the problems are filled with operators and coupling operators as well as
* boundary conditions are added.
*
* In the adaption loop the problems are solved the same order as they are added to the
* iterationInterface in this method. This order can be changed manually in the oneIteration
* method.
**/ **/
void initialize(AdaptInfo *adaptInfo) void initData()
{ {
for (size_t i = 0; i < chProb->getNumProblems(); i++)
addProblem(chProb->getProblem(i));
for (size_t i = 0; i < nsProb->getNumProblems(); i++)
addProblem(nsProb->getProblem(i));
addIterationInterface(chProb);
addIterationInterface(nsProb);
// setSolveProblem("ch",false);
addTimeInterface(chProb);
addTimeInterface(nsProb);
CouplingProblemStat::initialize(INIT_ALL);
dim = getMesh()->getDim();
fillCouplingOperators();
// fillOperators and fillBoundaryConditions for chProb and nsProb
nsProb->initTimeInterface();
chProb->initTimeInterface();
y = new DOFVector<double>(chProb->getFeSpace(0), "volume"); y = new DOFVector<double>(chProb->getFeSpace(0), "volume");
yc1 = new DOFVector<double>(chProb->getFeSpace(0), "volume"); yc1 = new DOFVector<double>(chProb->getFeSpace(0), "volume");
yc2 = new DOFVector<double>(chProb->getFeSpace(0), "volume"); yc2 = new DOFVector<double>(chProb->getFeSpace(0), "volume");
output = fopen( "output_data", "w" ); output = fopen( "output_data", "w" );
fprintf(output, "time u_b center_old\n"); fprintf(output, "time u_b center_old\n");
fclose(output); fclose(output);
u_b=0; u_b = 0;
// fillCouplingBoundaryConditions();
}
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{ FUNCNAME("CouplingIterationInterface::oneIteration()");
Flag flag = CouplingIterationInterface::oneIteration(adaptInfo, toDo);
//chProb->solveInitialProblem(adaptInfo);
return flag; // meshFunction for k local refinement around the interface of the phasefield-DOFVector
}; refFunction = new PhaseFieldRefinement(chProb->getMesh());
if (chProb->getDoubleWellType() == 0) {
refinement = new RefinementLevelDOF(
chProb->getFeSpace(),
refFunction,
chProb->getSolution()->getDOFVector(1)); // phaseField-DOFVector in [0,1]
} else {
refinement = new RefinementLevelDOF(
chProb->getFeSpace(),
refFunction,
new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(1))); // phaseField-DOFVector in [-1,1]
}
}
/** /**
...@@ -166,18 +135,21 @@ public: ...@@ -166,18 +135,21 @@ public:
void fillCouplingOperators() void fillCouplingOperators()
{ FUNCNAME("CahnHilliardNavierStokes::fillCouplingOperators()"); { FUNCNAME("CahnHilliardNavierStokes::fillCouplingOperators()");
MSG("CahnHilliardNavierStokes::fillCouplingOperators()"); MSG("CahnHilliardNavierStokes::fillCouplingOperators()");
int degree = nsProb->getFeSpace(0)->getBasisFcts()->getDegree();
for (size_t i = 0; i < dow; i++) { for (size_t i = 0; i < super::dow; i++) {
// < nu * d_i(c) , theta > // < nu * d_i(c) , theta >
Operator *opNuGradC = new Operator(nsProb->getFeSpace(i), chProb->getFeSpace(0)); Operator *opNuGradC = new Operator(nsProb->getFeSpace(i), chProb->getFeSpace(0));
opNuGradC->addTerm(new VecAndPartialDerivative_ZOT(chProb->getSolution()->getDOFVector(0), chProb->getSolution()->getDOFVector(1), i)); opNuGradC->addTerm(new VecAndPartialDerivative_ZOT(chProb->getSolution()->getDOFVector(0),
chProb->getSolution()->getDOFVector(1), i));
nsProb->getProblem(0)->addVectorOperator(opNuGradC, i, &surfaceTension, &surfaceTension); nsProb->getProblem(0)->addVectorOperator(opNuGradC, i, &surfaceTension, &surfaceTension);
} }
for (size_t i = 0; i < dow; i++) { for (size_t i = 0; i < super::dow; i++) {
// < v * grad(c) , theta > // < v * grad(c) , theta >
Operator *opVGradC = new Operator(chProb->getFeSpace(0), chProb->getFeSpace(1)); Operator *opVGradC = new Operator(chProb->getFeSpace(0), chProb->getFeSpace(1));
opVGradC->addTerm(new VecAtQP_FOT(nsProb->getSolution()->getDOFVector(i), new ID(), i), GRD_PHI); opVGradC->addTerm(new VecAndPartialDerivative_FOT(nsProb->getSolution()->getDOFVector(i), i), GRD_PHI);
chProb->getProblem()->addMatrixOperator(opVGradC, 1, 1, chProb->getTau()); chProb->getProblem()->addMatrixOperator(opVGradC, 1, 1, chProb->getTau());
} }
/**/ /**/
...@@ -191,27 +163,12 @@ public: ...@@ -191,27 +163,12 @@ public:
**/ **/
void solveInitialProblem(AdaptInfo *adaptInfo) void solveInitialProblem(AdaptInfo *adaptInfo)
{ {
// meshFunction for klocal refinement around the interface of the phasefield-DOFVector
refFunction = new PhaseFieldRefinement(chProb->getMesh());
if (chProb->getDoubleWellType() == 0) {
refinement = new RefinementLevelDOF(
chProb->getFeSpace(),
refFunction,
chProb->getSolution()->getDOFVector(1)); // phaseField-DOFVector in [0,1]
} else {
refinement = new RefinementLevelDOF(
chProb->getFeSpace(),
refFunction,
new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(1))); // phaseField-DOFVector in [-1,1]
}
// initial refinement // initial refinement
refinement->refine(0); refinement->refine(0);
for (int i = 0; i < 3; i++) { for (int i = 0; i < 5; i++) {
chProb->solveInitialProblem(adaptInfo); // update phaseField-DOFVector chProb->solveInitialProblem(adaptInfo); // update phaseField-DOFVector
refinement->refine((i < 4 ? 4 : 10)); // do k steps of local refinement refinement->refine(5); // do k steps of local refinement
} }
// solve all initial problems of the problems added to the CouplingTimeInterface // solve all initial problems of the problems added to the CouplingTimeInterface
...@@ -226,16 +183,14 @@ public: ...@@ -226,16 +183,14 @@ public:
void closeTimestep(AdaptInfo *adaptInfo) void closeTimestep(AdaptInfo *adaptInfo)
{ {
CouplingTimeInterface::closeTimestep(adaptInfo); CouplingTimeInterface::closeTimestep(adaptInfo);
AMDiS::Component<1> yFun(y->getFeSpace()->getBasisFcts()->getDegree());
Y yFun;
y->interpol(&yFun); y->interpol(&yFun);
DOFVector<double>::Iterator phiIt(dynamic_cast<DOFIndexed<double>*>(const_cast<DOFVector<double>*>(chProb->getSolution()->getDOFVector(1))), USED_DOFS); DOFIterator<double> phiIt(chProb->getSolution()->getDOFVector(1), USED_DOFS);
DOFVector<double>::Iterator yIt(dynamic_cast<DOFIndexed<double>*>(y), USED_DOFS); DOFIterator<double> yIt(y, USED_DOFS);
DOFVector<double>::Iterator yc1It(dynamic_cast<DOFIndexed<double>*>(yc1), USED_DOFS); DOFIterator<double> yc1It(yc1, USED_DOFS);
DOFVector<double>::Iterator yc2It(dynamic_cast<DOFIndexed<double>*>(yc2), USED_DOFS); DOFIterator<double> yc2It(yc2, USED_DOFS);
for (phiIt.reset(), yIt.reset(),yc1It.reset(),yc2It.reset(); for (phiIt.reset(), yIt.reset(),yc1It.reset(),yc2It.reset();
!phiIt.end(); !phiIt.end();
...@@ -244,25 +199,18 @@ public: ...@@ -244,25 +199,18 @@ public:
*yc1It = *yc2It * *yIt; *yc1It = *yc2It * *yIt;
} }
if (adaptInfo->getTimestepNumber() <= 1) center_old = (yc1->Int())/(yc2->Int()); if (adaptInfo->getTimestepNumber() <= 1)
u_b = - (center_old - (yc1->Int())/(yc2->Int()) ) / adaptInfo->getTimestep(); center_old = (yc1->Int()) / (yc2->Int());
center_old = (yc1->Int())/(yc2->Int()); u_b = - (center_old - (yc1->Int()) / (yc2->Int()) ) / adaptInfo->getTimestep();
center_old = (yc1->Int()) / (yc2->Int());
// append values to file
output = fopen( "output_data", "a" ); output = fopen( "output_data", "a" );
fprintf(output, "%f, %f, %f\n", fprintf(output, "%f, %f, %f\n", adaptInfo->getTime(), u_b, center_old);
adaptInfo->getTime(),
u_b,
center_old);
fclose(output); fclose(output);
//chProb->solveInitialProblem(adaptInfo);
refinement->refine(2); refinement->refine(2);
} }
protected: protected:
...@@ -272,9 +220,6 @@ protected: ...@@ -272,9 +220,6 @@ protected:
PhaseFieldRefinement* refFunction; PhaseFieldRefinement* refFunction;
RefinementLevelDOF *refinement; RefinementLevelDOF *refinement;
unsigned dim; // dimension of the mesh
unsigned dow; // dimension of the world
double sigma; // coupling parameter to calculate the surface tension double sigma; // coupling parameter to calculate the surface tension
double surfaceTension;// := sigma/epsilon double surfaceTension;// := sigma/epsilon
double u_b, center_old; double u_b, center_old;
......
...@@ -115,11 +115,11 @@ void CahnHilliard_::solveInitialProblem(AdaptInfo *adaptInfo) ...@@ -115,11 +115,11 @@ void CahnHilliard_::solveInitialProblem(AdaptInfo *adaptInfo)
/// create phase-field from signed-dist-function /// create phase-field from signed-dist-function
if (doubleWell == 0) { if (doubleWell == 0) {
transformDOF(prob->getSolution()->getDOFVector(1), forEachDOF(prob->getSolution()->getDOFVector(1),
new SignedDistToPhaseField(initialEps)); new SignedDistToPhaseField(initialEps, 1.0/(2.0*sqrt(2.0))));
} else { } else {
transformDOF(prob->getSolution()->getDOFVector(1), forEachDOF(prob->getSolution()->getDOFVector(1),
new SignedDistToCh(initialEps)); new SignedDistToCh(initialEps, 1.0/sqrt(2.0)));
} }
} }
} }
...@@ -130,6 +130,7 @@ void CahnHilliard_::fillOperators() ...@@ -130,6 +130,7 @@ void CahnHilliard_::fillOperators()
MSG("CahnHilliard_::fillOperators()\n");