From 5be7bb7e64cb8cda1474a6d908be3c3da4d7f288 Mon Sep 17 00:00:00 2001 From: Simon Praetorius <simon.praetorius@tu-dresden.de> Date: Tue, 27 Mar 2012 17:46:45 +0000 Subject: [PATCH] some new domes --- demo/src/chns/CahnHilliardNavierStokes.h | 167 +++++++++++++++++++++++ demo/src/chns/drivenCavity.cc | 84 ++++++++++++ 2 files changed, 251 insertions(+) create mode 100644 demo/src/chns/CahnHilliardNavierStokes.h create mode 100644 demo/src/chns/drivenCavity.cc diff --git a/demo/src/chns/CahnHilliardNavierStokes.h b/demo/src/chns/CahnHilliardNavierStokes.h new file mode 100644 index 00000000..39d7c407 --- /dev/null +++ b/demo/src/chns/CahnHilliardNavierStokes.h @@ -0,0 +1,167 @@ +/** \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 diff --git a/demo/src/chns/drivenCavity.cc b/demo/src/chns/drivenCavity.cc new file mode 100644 index 00000000..8a70118f --- /dev/null +++ b/demo/src/chns/drivenCavity.cc @@ -0,0 +1,84 @@ +#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; +}; -- GitLab