Commit 5be7bb7e authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

some new domes

parent 7f3072a4
/** \file CahnHilliardNavierStokes.h */
#ifndef CAHN_HILLIARD_NAVIER_STOKES_H
#define CAHN_HILLIARD_NAVIER_STOKES_H
// coupling structures
#include "CouplingIterationInterface.h"
#include "CouplingTimeInterface.h"
#include "CouplingProblemStat.h"
// structures for local refinement
#include "Refinement.h"
#include "MeshFunction_Level.h"
#include "POperators.h"
#include "CouplingOperators.h"
using namespace AMDiS;
/**
* \ingroup Problem
*
* \brief
*/
template<typename CH_Type=CahnHilliard, typename NS_Type=NavierStokes_TaylorHood>
class CahnHilliardNavierStokes : public CouplingIterationInterface,
public CouplingTimeInterface,
public CouplingProblemStat
{
public:
CahnHilliardNavierStokes(std::string name_, CH_Type *chProb_, NS_Type *nsProb_)
: CouplingProblemStat(name_),
chProb(chProb_),
nsProb(nsProb_),
refFunction(NULL),
refinement(NULL),
sigma(0.0),
surfaceTension(0.0)
{
dow = Global::getGeo(WORLD);
Parameters::get(name + "->sigma", sigma);
surfaceTension = sigma/chProb->getEpsilon();
}
~CahnHilliardNavierStokes() {
if (refFunction != NULL)
delete refFunction;
if (refinement != NULL)
delete refinement;
}
void initialize(AdaptInfo *adaptInfo)
{
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);
addTimeInterface(chProb);
addTimeInterface(nsProb);
CouplingProblemStat::initialize(INIT_ALL);
dim = getMesh()->getDim();
fillCouplingOperators();
// fillOperators and fillBoundaryConditions for chProb and nsProb
nsProb->initTimeInterface();
chProb->initTimeInterface();
fillCouplingBoundaryConditions();
}
void fillCouplingOperators()
{ FUNCNAME("CahnHilliardNavierStokes::fillCouplingOperators()");
MSG("CahnHilliardNavierStokes::fillCouplingOperators()");
for (size_t i = 0; i < dow; i++) {
// < nu * d_i(c) , theta >
Operator *opNuGradC = new Operator(nsProb->getFeSpace(i), chProb->getFeSpace(0));
opNuGradC->addTerm(new VecAndPartialDerivative_ZOT(chProb->getSolution()->getDOFVector(1), chProb->getSolution()->getDOFVector(0), i));
nsProb->getProblem(0)->addVectorOperator(opNuGradC, i, &surfaceTension, &surfaceTension);
//
// // stabilizing term
// Operator *stabilBeltrami = new Operator(nsProb->getFeSpace(i), nsProb->getFeSpace(i));
// stabilBeltrami->addTerm(new MatrixGradient_SOT(
// chProb->getSolution()->getDOFVector(0), new ProjectionMatrix(surfaceTension), new DivFct())); // factor = sigma ?
// nsProb->getProblem(0)->addMatrixOperator(stabilBeltrami, i, i, nsProb->getTau(), nsProb->getTau());
}
// < v * grad(c) , theta >
Operator *opVGradC = new Operator(chProb->getFeSpace(0), chProb->getFeSpace(0));
opVGradC->addTerm(new WorldVector_FOT(nsProb->getVelocity(), -1.0), GRD_PSI);
chProb->getProblem()->addMatrixOperator(opVGradC, 0, 0);
}
void fillCouplingBoundaryConditions()
{ FUNCNAME("CahnHilliardNavierStokes::fillCouplingBoundaryConditions()");
}
/// Solves the initial problem.
virtual void solveInitialProblem(AdaptInfo *adaptInfo)
{
refFunction= new PhaseFieldRefinement(chProb->getMesh());
refinement= new RefinementLevelDOF(
chProb->getFeSpace(),
refFunction,
new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(0))); // phaseField-DOFVector
// set initial values
refinement->refine(0);
for (int i= 0; i < 3; ++i) {
chProb->solveInitialProblem(adaptInfo); // initial phaseField
refinement->refine((i < 4 ? 4 : 10));
}
CouplingTimeInterface::solveInitialProblem(adaptInfo);
}
/// Called at the beginning of each timestep
virtual void initTimestep(AdaptInfo *adaptInfo)
{
CouplingTimeInterface::initTimestep(adaptInfo);
}
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
{
CouplingIterationInterface::oneIteration(adaptInfo, toDo);
}
/// Called at the end of each timestep.
virtual void closeTimestep(AdaptInfo *adaptInfo)
{
CouplingTimeInterface::closeTimestep(adaptInfo);
refinement->refine(2);
}
protected:
CH_Type *chProb;
NS_Type *nsProb;
PhaseFieldRefinement* refFunction;
RefinementLevelDOF *refinement;
unsigned dim;
unsigned dow;
double sigma;
double surfaceTension;
};
#endif // CAHN_HILLIARD_NAVIER_STOKES_H
#include "AMDiS.h"
#include "CahnHilliard.h"
#include "NavierStokes_TaylorHood.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;
struct DrivenCavityBC : AbstractFunction<double, WorldVector<double> >
{
double operator()(const WorldVector<double> &x) const {
return std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
}
};
class CH_DrivenCavity : public CahnHilliard
{
public:
CH_DrivenCavity(std::string name_) : CahnHilliard(name_) {}
protected:
void fillBoundaryConditions()
{ FUNCNAME("CH_DrivenCavity::fillBoundaryConditions()");
// homogeneouse neumann conditions
}
};
class NS_DrivenCavity : public NavierStokes_TaylorHood
{
public:
NS_DrivenCavity(std::string name_) : NavierStokes_TaylorHood(name_) {}
protected:
void fillBoundaryConditions()
{ FUNCNAME("NS_DrivenCavity::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
for (size_t i = 0; i < dow; i++)
getProblem(0)->addDirichletBC(1, i, i, zeroDOF);
/// at upper wall: prescribed velocity
getProblem(0)->addDirichletBC(2, 0, 0, new DrivenCavityBC);
getProblem(0)->addDirichletBC(2, 1, 1, zeroDOF);
}
};
int main(int argc, char** argv)
{ FUNCNAME("main");
AMDiS::init(argc, argv);
CH_DrivenCavity chProb("ch");
NS_DrivenCavity nsProb("ns");
CahnHilliardNavierStokes<CH_DrivenCavity,NS_DrivenCavity> mainProb("main", &chProb, &nsProb);
// Adapt-Infos
AdaptInfo adaptInfo("adapt", mainProb.getNumComponents());
mainProb.initialize(&adaptInfo);
// adaption loop - solve ch-prob and ns-prob
AdaptInstationary adaptInstat("adapt", mainProb, adaptInfo, mainProb, 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;
};
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment