diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index b032d20434f0f1075f6fe8136e8897671d8dc5f0..104a8911a3ea863e5b456da1fcdcaf8e1ce706ab 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -29,24 +29,121 @@
//_________________________________________________________________________________________________________________________
-template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+template<typename Grid>
class chns_base
{
public:
typedef typename Grid::LeafGridView GV;
const int dim = GV::dimension;
- GV gv;
typedef double RF;
chns_base(Grid& grid, Dune::ParameterTree& ptree)
{
- GV gv = grid.leafGridView();
}
};
-template<typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+class TimestepManager
+{
+public:
+ double T;
+ double dt;
+ double t;
+
+ double dtMin;
+ double dtMax;
+ double dtIncreaseFactor;
+ double dtDecreaseFactor;
+
+ double lastsave;
+ double saveintervall;
+
+ int verb;
+
+ TimestepManager(const Dune::ParameterTree & param, int verbosity)
+ {
+ verb = verbosity;
+ T=param.get("tMax",(double)1.0);
+ dt=param.get("dt",(double)0.1);
+
+ dtMin = param.get("dtMin",(double)dt);
+ dtMax = param.get("dtMax",(double)dt);
+ dtIncreaseFactor = param.get("dtIncreaseFactor",(double)1.2);
+ dtDecreaseFactor = param.get("dtDecreaseFactor",(double)0.5);
+
+ saveintervall = param.get("temporal.saveintervall",(double)0.0);
+
+ t=param.get("tIni",(double)0.0);
+ lastsave = t - saveintervall;
+ }
+
+ // RETURN VALUE: -1 -> ABORT, 0-> REDO TIMESTEP, 1-> OK
+ template<typename NS_NEWTON, typename CH_NEWTON, typename VTKWRITER>
+ int doTimestep(NS_NEWTON& ns_newton, CH_NEWTON& ch_newton, VTKWRITER& vtkwriter)
+ {
+ try
+ {
+ if (verb > 1) std::cout << "= Begin newtonapply CH:" << std::endl;
+ ch_newton.apply();
+
+ if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
+ ns_newton.apply();
+ }
+ catch(...)
+ {
+ if(dt <= dtMin)
+ {
+ std::cout << "===============================" << std::endl;
+ if(dtMin==dtMax)
+ std::cout << "Problem with Newton Solver, aborting :(" << std::endl;
+ else
+ std::cout << "Timestep too small, aborting :(" << std::endl;
+ std::cout << "===============================" << std::endl;
+ return -1;
+ }
+ dt = std::max(dt*dtDecreaseFactor,dtMin);
+ if (verb > 0) std::cout << "Reducing Timestep to "<< dt << std::endl;
+
+ return 0;
+ }
+
+ t += dt;
+ if(dt < dtMax)
+ {
+ dt = std::min(dt*dtIncreaseFactor, dtMax);
+ if (verb > 0) std::cout << "New Timestep "<< dt << std::endl;
+ }
+
+ if( t>lastsave+saveintervall )
+ {
+ lastsave = t;
+ vtkwriter.write(t,Dune::VTK::ascii);
+ }
+ return 1;
+ }
+
+
+ bool isFinished()
+ {
+ return (t>=T);
+ }
+};
+
+template<typename GFS, typename NS_BDRY_FN, typename CH_BDRY_FN, typename NS_V, typename CH_V>
+void interpolateToBdry(const GFS& gfs, const NS_BDRY_FN& nsBdryFn, const CH_BDRY_FN& chBdryFn, NS_V& nsV, CH_V& chV)
+{
+ NS_V nsInterp(gfs.ns);
+ CH_V chInterp(gfs.ch);
+ Dune::PDELab::interpolate(nsBdryFn, gfs.ns, nsInterp);
+ Dune::PDELab::interpolate(chBdryFn, gfs.ch, chInterp);
+ Dune::PDELab::copy_constrained_dofs(gfs.nsCC, nsInterp, nsV);
+ Dune::PDELab::copy_constrained_dofs(gfs.chCC, chInterp, chV);
+}
+
+//_________________________________________________________________________________________________________________________
+
+template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
class GFS_2f0s
{
public:
@@ -56,8 +153,10 @@ public:
typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
- typedef Dune::PDELab::VectorGridFunctionSpace
- <ES,V_FEM,ES::GridView::dimension,VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+ ConstraintsAssembler conp, conphi, conmu;
+
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
+ VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
V_GFS vGfs;
typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
@@ -69,399 +168,206 @@ public:
typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> MU_GFS;
MU_GFS muGfs;
- typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS_GFS;
- NS_GFS nsGfs;
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS;
+ NS ns;
- typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS> CH_GFS;
- CH_GFS chGfs;
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS> CH;
+ CH ch;
- GFS_2f0s(const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem)
- {
- ConstraintsAssembler conp, conphi, conmu;
+ typedef typename NS::template ConstraintsContainer<RF>::Type NS_CC;
+ NS_CC nsCC;
- vGfs(es,vFem);
- vGfs.name("velocity");
+ typedef typename CH::template ConstraintsContainer<RF>::Type CH_CC;
+ CH_CC chCC;
- pGfs(es,pFem, conp);
+ GFS_2f0s(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem) :
+ conp(),conphi(),conmu(),
+ vGfs(es,vFem),
+ pGfs(es,pFem, conp),
+ phiGfs(es,phiFem, conphi),
+ muGfs(es,phiFem, conmu),
+ ns(vGfs, pGfs),
+ ch(phiGfs, muGfs),
+ nsCC(),
+ chCC()
+ {
+ vGfs.name("velocity");
pGfs.name("pressure");
-
- phiGfs(es,phiFem, conphi);
phiGfs.name("phasefield");
-
- muGfs(es,phiFem, conmu);
muGfs.name("potential");
+ }
+
+ template<typename Bdry>
+ void updateConstraintsContainer(const Bdry& bdry)
+ {
+ Dune::PDELab::constraints(bdry.bdryNS,ns,nsCC, 0);
+ Dune::PDELab::constraints(bdry.bdryCH,ch,chCC, 0);
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+
+template<typename GFS, typename NS_V, typename CH_V>
+class DGF_2f0s
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VSubGFS;
+ VSubGFS vSubGfs;
+ typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,NS_V> VDGF;
+ VDGF vDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<1> > PSubGFS;
+ PSubGFS pSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,NS_V> PDGF;
+ PDGF pDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<0> > PhiSubGFS;
+ PhiSubGFS phiSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PhiSubGFS,CH_V> PhiDGF;
+ PhiDGF phiDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<1> > MuSubGFS;
+ MuSubGFS muSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<MuSubGFS,CH_V> MuDGF;
+ MuDGF muDgf;
+
- nsGfs(vGfs, pGfs);
- chGfs(phiGfs, muGfs);
+ DGF_2f0s(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
+ vSubGfs(gfs.ns),
+ vDgf(vSubGfs,nsV),
+ pSubGfs(gfs.ns),
+ pDgf(pSubGfs,nsV),
+ phiSubGfs(gfs.ch),
+ phiDgf(phiSubGfs,chV),
+ muSubGfs(gfs.ch),
+ muDgf(muSubGfs,chV)
+ {
+ }
+
+ template<typename VTKWRITER>
+ void addToVTKWriter(VTKWRITER& vtkwriter)
+ {
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vDgf,"v"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
}
};
-// void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
-// {
-//
-//
-// //std::bitset<2> periodic_directions(std::string("01"));
-// //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
-//
-// // dimension and important types
-// // type for computations
-//
-// RF T = ptree.get("temporal.time",(RF)1.0);
-//
-//
-// RF initialvmax = ptree.get("initial.vmax",(RF)1);
-// RF vmaxtime = ptree.get("initial.vmaxtime",(RF)1);
-//
-// RF time = 0.0;
-// RF lastdatasave = 0.0;
-// RF datasaveintervall = ptree.get("temporal.saveintervall",(RF)1.0);
-// RF dt_min = 1e-8;
-// bool savealways = ptree.get("output.savealways",(RF)1);
-//
-// Problem<RF> problem(ptree.get("temporal.tau",(RF)0.1),ptree);
-// problem.setTime(time);
-//
-//
-// //_________________________________________________________________________________________________________________________
-// //Initial Values
-// //typedef ZeroScalarFunction<GV,RF> InitialPressure;
-// //InitialPressure initialPressure(gv);
-// typedef TUPressure<GV,RF,1> InitialPressure;
-// InitialPressure initialPressure(gv);
-//
-// typedef IniVelocity<GV,RF,2> InitialVelocity;
-// InitialVelocity initialVelocity(gv,initialvmax,vmaxtime);
-// typedef Dune::PDELab::CompositeGridFunction<InitialVelocity,InitialPressure> InitialNS;
-// InitialNS initialNS(initialVelocity,initialPressure);
-//
-// typedef PhiZero<GV,RF,1> InitialPhi;
-// InitialPhi initialPhi(gv,0.5,0.15,problem.eps);
-// typedef ThreephasePhiZero<GV,RF,1> InitialThreephasePhi;
-// InitialThreephasePhi initialPhiOne(gv,0, 0.5,0.15,problem.eps);
-// InitialThreephasePhi initialPhiTwo(gv,1, 0.5,0.15,problem.eps);
-// typedef ZeroScalarFunction<GV,RF> InitialMu;
-// InitialMu initialMu(gv);
-//
-// typedef Dune::PDELab::CompositeGridFunction<InitialPhi,InitialMu> InitialCH;
-// InitialCH initialCH(initialPhi,initialMu);
-//
-//
-// //_________________________________________________________________________________________________________________________
-// // Make grid function spaces
-//
-//
-//
-//
-// //_________________________________________________________________________________________________________________________
-// //Constraints
-// //typedef Dune::PDELab::EmptyTransformation CC;
-// typedef typename NS_GFS::template ConstraintsContainer<RF>::Type NSC;
-// NSC nscg;
-// nscg.clear();
-//
-// typedef typename CH_GFS::template ConstraintsContainer<RF>::Type CHC;
-// CHC chcg;
-//
-// // create Taylor-Hood constraints (no-slip, no penetration)
-// typedef Dune::PDELab::DirichletConstraintsParameters AlwaysDirichlet;
-// typedef Dune::PDELab::NoDirichletConstraintsParameters NeverDirichlet;
-// typedef Dune::PDELab::PowerConstraintsParameters<AlwaysDirichlet,dim> VelocityConstraints;
-// typedef Dune::PDELab::CompositeConstraintsParameters<VelocityConstraints, SinglePointPressureBC> NSConstraints;
-//
-// AlwaysDirichlet scalarvelocity_constraints;
-// SinglePointPressureBC pressure_constraints;
-// VelocityConstraints velocity_constraints(scalarvelocity_constraints);
-// NSConstraints nsConstraints(velocity_constraints,pressure_constraints);
-//
-// NeverDirichlet pf_constraints;
-// NeverDirichlet mu_constraints;
-//
-// typedef Dune::PDELab::CompositeConstraintsParameters<NeverDirichlet, NeverDirichlet> CHChonstraints;
-// CHChonstraints chConstraints(pf_constraints,mu_constraints);
-//
-//
-//
-// Dune::PDELab::constraints(nsConstraints,nsGfs,nscg, 0);
-// Dune::PDELab::constraints(chConstraints,chGfs,chcg, 0);
-//
-//
-// //_________________________________________________________________________________________________________________________
-// //Coefficient vectors
-//
-// using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
-// using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
-// CH_V chV(chGfs);
-// NS_V nsV(nsGfs);
-// chV=1.1245;
-// nsV=2;
-// CH_V chVOld = chV;
-// NS_V nsVOld = nsV;
-// CH_V chInterpolant = chV;
-// NS_V nsInterpolant = nsV;
-//
-// //_________________________________________________________________________________________________________________________
-// // LocalOperators
-//
-// typedef CHNS_FEM_NavierStokes<Problem<RF>, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
-// NS_LOP nsLop(problem, chGfs, chV, chVOld, nsGfs, nsVOld);
-//
-// typedef CHNS_FEM_CahnHilliard<Problem<RF>, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
-// CH_LOP chLop(problem, chGfs, chVOld, nsGfs, nsV);
-//
-// typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
-// MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
-//
-//
-//
-//
-//
-// //_________________________________________________________________________________________________________________________
-// // Interpolate Initial Data
-// //Dune::PDELab::interpolate(initialFn,nsGfs,nsVOld);
-//
-// Dune::PDELab::interpolate(initialCH,chGfs,chV);
-// Dune::PDELab::interpolate(initialNS,nsGfs,nsV);
-//
-//
-//
-// //_________________________________________________________________________________________________________________________
-// // Construct discrete Grid Functions
-// typedef typename Dune::PDELab::GridFunctionSubSpace<NS_GFS,Dune::TypeTree::TreePath<0> > VelocitySubGFS;
-// VelocitySubGFS velocitySubGfs(nsGfs);
-// typedef typename Dune::PDELab::GridFunctionSubSpace<NS_GFS,Dune::TypeTree::TreePath<1> > PressureSubGFS;
-// PressureSubGFS pressureSubGfs(nsGfs);
-// typedef Dune::PDELab::VectorDiscreteGridFunction<VelocitySubGFS,NS_V> VDGF;
-// VDGF vDgf(velocitySubGfs,nsV);
-// typedef Dune::PDELab::DiscreteGridFunction<PressureSubGFS,NS_V> PDGF;
-// PDGF pDgf(pressureSubGfs,nsV);
-//
-//
-// typedef typename Dune::PDELab::GridFunctionSubSpace<CH_GFS,Dune::TypeTree::TreePath<0> > PhasefieldSubGFS;
-// PhasefieldSubGFS phasefieldSubGfs(chGfs);
-//
-// typedef typename Dune::PDELab::GridFunctionSubSpace<CH_GFS,Dune::TypeTree::TreePath<1> > PotentialSubGFS;
-// PotentialSubGFS potentialSubGfs(chGfs);
-//
-// typedef Dune::PDELab::DiscreteGridFunction<PotentialSubGFS,CH_V> MuDGF;
-// MuDGF muDgf(potentialSubGfs,chV);
-//
-//
-// typedef Dune::PDELab::DiscreteGridFunction<PhasefieldSubGFS,CH_V> PhiDGF;
-// PhiDGF phiDgf(phasefieldSubGfs,chV);
-//
-// //_________________________________________________________________________________________________________________________
-// // prepare VTK writer
-// std::string filename=ptree.get("output.filename","output");
-// struct stat st;
-// if( stat( filename.c_str(), &st ) != 0 )
-// {
-// int stat = 0;
-// stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
-// if( stat != 0 && stat != -1)
-// std::cout << "Error: Cannot create directory "
-// << filename << std::endl;
-// }
-// auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
-// Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
-// vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vDgf,"v"));
-// vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
-// vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
-// vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
-//
-//
-// //_________________________________________________________________________________________________________________________
-// int i=0;
-// for(auto it = chV.begin(); it!=chV.end(); ++it)
-// {
-// ++i;
-// if(i<100) std::cout<< *it <<std::endl;
-// }
-// std::cout << i << std::endl;
-// i=0;
-// for(auto it = nsV.begin(); it!=nsV.end(); ++it)
-// {
-// ++i;
-// if(i<100) std::cout<< *it <<std::endl;
-// }
-// std::cout << i << std::endl;
-//
-// Dune::PDELab::set_constrained_dofs(nscg,0.0,nsV);
-// //Dune::PDELab::set_constrained_dofs(chcg,0.0,chV);
-//
-// vtkwriter.write(time,Dune::VTK::ascii);
-//
-//
-// //_________________________________________________________________________________________________________________________
-// //int i_it=0;
-// chVOld = chV;
-// nsVOld = nsV;
-// while (time < T - dt_min*0.5)
-// {
-// problem.setTime(time+problem.dt);
-//
-// // Grid Operators
-// typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,NSC,NSC> NS_GO;
-// NS_GO nsGo(nsGfs,nscg,nsGfs,nscg,nsLop,mbe);
-//
-// typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,CHC,CHC> CH_GO;
-// CH_GO chGo(chGfs,chcg,chGfs,chcg,chLop,mbe);
-//
-// // Linear solver
-// //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 CH_LS; CH_LS ls_ch(200,5000,1);
-// //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_SSORk<CH_GO> CH_LS; CH_LS ls_ch(chGo,5000,3,1);
-//
-// typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
-// <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,3,1);
-//
-//
-// //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 NS_LS; NS_LS ls_ns(200,5000,1);
-// //typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_ILU0 NS_LS; NS_LS ls_ns(5000,1);
-//
-// //AT LEAST WORKS:
-// //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_NOPREC<NS_GFS> NS_LS; NS_LS ls_ns(nsGfs,5000,1);
-//
-// //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_Jacobi<NS_GFS> NS_LS; NS_LS ls_ns(nsGfs,5000,1);
-// typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
-// <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
-// //typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_Jac LS;
-// //typedef Dune::PDELab::ISTLBackend_SEQ_SuperLU LS;
-// //typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LS;
-//
-// // Solve (possibly) nonlinear problem
-// typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
-// PDESolver1 ns_newton(nsGo,nsV,ls_ns);
-// ns_newton.setReassembleThreshold(0.0);
-// ns_newton.setVerbosityLevel(2);
-// ns_newton.setMaxIterations(50);
-// ns_newton.setLineSearchMaxIterations(30);
-//
-// ns_newton.setAbsoluteLimit(1e-9);
-// ns_newton.setReduction(1e-7);
-//
-// typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
-// PDESolver2 ch_newton(chGo,chV,ls_ch);
-// ch_newton.setReassembleThreshold(0.0);
-// ch_newton.setVerbosityLevel(1);
-// ch_newton.setMaxIterations(50);
-// ch_newton.setLineSearchMaxIterations(50);
-//
-//
-//
-// //_________________________________________________________________________________________________________________________
-// // do time step
-// initialNS.setTime(time);
-// Dune::PDELab::interpolate(initialNS,nsGfs,nsInterpolant);
-// Dune::PDELab::interpolate(initialCH,chGfs,chInterpolant);
-// Dune::PDELab::copy_constrained_dofs(nscg, nsInterpolant, nsV);
-// Dune::PDELab::copy_constrained_dofs(chcg, chInterpolant, chV );
-// bool EverythingOk = true;
-// try
-// {
-// if (helper.rank()==0)
-// std::cout << "=== Begin newtonapply CH:" << std::endl;
-// ch_newton.apply();
-// //vtkwriter.write(time+5e-6,Dune::VTK::ascii);
-// if (helper.rank()==0)
-// std::cout << "=== Begin newtonapply NS:" << std::endl;
-// ns_newton.apply();
-// }
-// catch(...)
-// {
-// //vtkwriter.write(time+6e-6,Dune::VTK::ascii);
-// if (helper.rank()==0)
-// std::cout << "Reducing Timestep to "<< problem.dt/2 << std::endl;
-// problem.dt = problem.dt/2;
-// EverythingOk = false;
-// //time+=1e-5;
-// //NS_V nsRes(nsGfs);
-// //nsGo.residual(nsV,nsRes);
-// //nsV = nsRes;
-// //vtkwriter.write(time,Dune::VTK::ascii);
-// nsV = nsVOld;
-// chV = chVOld;
-// if(problem.dt<ptree.get("temporal.taumin",(RF)0.0001))
-// {
-// if (helper.rank()==0)
-// std::cout << "Timestep too small, aborting :(" << std::endl;
-// return;
-// }
-// }
-//
-// if(EverythingOk)
-// {
-// time += problem.dt;
-// problem.dt = fmin(1.03*problem.dt, ptree.get("temporal.taumax",(RF)1000) );
-// if (helper.rank()==0)
-// std::cout << "New Timestep "<< problem.dt << std::endl;
-// if( (time>lastdatasave+datasaveintervall) || savealways)
-// {
-// lastdatasave = time;
-// vtkwriter.write(time,Dune::VTK::ascii);
-// }
-// chVOld = chV;
-// nsVOld = nsV;
-// }
-//
-// //_________________________________________________________________________________________________________________________
-// /*
-// //Estimator GFS
-// typedef Dune::PDELab::P0LocalFiniteElementMap<RF,RF,dim> P0FEM;
-// P0FEM p0fem(Dune::GeometryTypes::simplex(dim));
-// typedef Dune::PDELab::NoConstraints NCON;
-// typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM,NCON,VectorBackend> P0GFS;
-// P0GFS p0Gfs(gv,p0fem);
-// //Estimator vector
-// using Z0 = Dune::PDELab::Backend::Vector<P0GFS,RF>;
-// Z0 z0(p0Gfs,0.0);
-// typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
-// Z0DGF z0dgf(p0Gfs,z0);
-// //Estimator Operator
-// typedef ACFS_FEM_AC_Estimator<Problem<RF>,CH_GFS> EST_LOP;
-// EST_LOP estlop(problem);
-// typedef Dune::PDELab::EmptyTransformation NCC;
-// typedef Dune::PDELab::GridOperator<CH_GFS,P0GFS,EST_LOP,MBE,RF,RF,RF,NCC,NCC> ESTGO;
-// ESTGO estgo(chGfs,p0Gfs,estlop,mbe);
-// estgo.residual(chV,z0);
-// // mark elements for refinement
-// std::cout << "mark elements" << std::endl;
-// RF eta_refine = ptree.get("fem.etarefine",(double)0.15);
-// RF eta_coarsen = ptree.get("fem.etarecoarsen",(double)0.0);
-// int Rmaxlevel = ptree.get("fem.maxlevel",(double)3);
-// int Rminlevel = ptree.get("fem.minlevel",(double)2);
-// Dune::PDELab::mark_grid(grid,z0,eta_refine,eta_coarsen, Rminlevel,Rmaxlevel,2);
-//
-// //znew = z;
-// //znew = 1.0;
-// std::cout << "Size Before:" << chGfs.globalSize() << std::endl;
-//
-// // do refinement
-// std::cout << "adapt grid and solution" << std::endl;
-// auto chTransfer = transferSolutions(chGfs,2,chV,chVOld);
-// auto nsTransfer = transferSolutions(nsGfs,2,nsV,nsVOld);
-// Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
-// //Dune::PDELab::adapt_grid(grid, chGfs, chV, 2*2);
-// std::cout << "adapt grid and solution ended" << std::endl;
-// //znew = z;
-// //znew = 1.0;
-// std::cout << "Size After:" << chGfs.globalSize() << std::endl;
-// // recompute constraints
-// Dune::PDELab::constraints(nsConstraints,nsGfs,nscg, 0);
-// Dune::PDELab::constraints(chConstraints,chGfs,chcg, 0);
-// Dune::PDELab::interpolate(initialNS,nsGfs,nsVOld);
-// Dune::PDELab::interpolate(initialCH,chGfs,chVOld);
-// if(i_it>5)
-// {
-// Dune::PDELab::copy_constrained_dofs(nscg, nsV , nsVOld );
-// Dune::PDELab::copy_constrained_dofs(chcg, chV , chVOld );
-// }
-// else
-// {
-// chV=chVOld;
-// nsV=nsVOld;
-// }
-// i_it++;
-//
-// vtkwriter.write(time+1e-6,Dune::VTK::ascii);*/
-// }
-// }
-// }
+//_________________________________________________________________________________________________________________________
+
+template<typename IniV, typename IniP, typename IniPhi, typename IniMu>
+class Ini_2f0s
+{
+public:
+ IniV iniV;
+ IniP iniP;
+ IniPhi iniPhi;
+ IniMu iniMu;
+ typedef Dune::PDELab::CompositeGridFunction<IniV,IniP> IniNS;
+ typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu> IniCH;
+ IniNS iniNS;
+ IniCH iniCH;
+
+ template<typename GV, typename Params>
+ Ini_2f0s(const GV& gv, const Params& params) :
+ iniV(gv, params),
+ iniP(gv, params),
+ iniPhi(gv, params),
+ iniMu(gv, params),
+ iniNS(iniV,iniP),
+ iniCH(iniPhi,iniMu)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<const int dim, typename BdryV_Scalar, typename BdryP, typename BdryPhi, typename BdryMu>
+class Bdry_2f0s
+{
+public:
+ BdryV_Scalar bdryV_Scalar;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ Bdry_2f0s() :
+ bdryV_Scalar(),
+ bdryP(),
+ bdryPhi(),
+ bdryMu(),
+ bdryV(bdryV_Scalar),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename GV, typename RF, typename EST_LOP, typename CH_GFS>
+class RefinementIndicator
+{
+public:
+ //Estimator GFS
+ typedef Dune::PDELab::P0LocalFiniteElementMap<RF,RF,GV::dimension> P0FEM;
+ P0FEM p0fem;
+ typedef Dune::PDELab::NoConstraints NCON;
+ typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM,NCON,Dune::PDELab::ISTL::VectorBackend<> > P0GFS;
+ P0GFS p0Gfs;
+ //Estimator vector
+ using Z0 = Dune::PDELab::Backend::Vector<P0GFS,RF>;
+ Z0 z0;
+
+ //typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
+ //Z0DGF z0dgf;
+
+ //Estimator Operator
+ EST_LOP estlop;
+ typedef Dune::PDELab::EmptyTransformation NCC;
+ typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+ typedef Dune::PDELab::GridOperator<CH_GFS,P0GFS,EST_LOP,MBE,RF,RF,RF,NCC,NCC> ESTGO;
+ ESTGO estgo;
+
+ double eta_refine;
+ double eta_coarsen;
+ int Rmaxlevel;
+ int Rminlevel;
+
+ RefinementIndicator(const GV& gv, const CH_GFS chGfs, const Dune::ParameterTree & param) :
+ p0fem(Dune::GeometryTypes::simplex(GV::dimension)),
+ p0Gfs(gv,p0fem),
+ z0(p0Gfs,0.0),
+ //z0dgf(p0Gfs,z0),
+ estlop(),
+ estgo(chGfs,p0Gfs,estlop,MBE(1))
+ {
+ eta_refine = param.get("etaRefine",(double)0.15);
+ eta_coarsen = param.get("etaCoarsen",(double)0.05);
+ Rmaxlevel = param.get("maxRefineLevel",(double)2);
+ Rminlevel = param.get("minRefineLevel",(double)0);
+ }
+
+ template<typename CH_V>
+ const Z0& calculateIndicator(const CH_V& chV)
+ {
+ estgo.residual(chV,z0);
+ return z0;
+ }
+
+ template<typename Grid>
+ void markGrid(Grid &grid, int verbosity=0)
+ {
+ if(verbosity>0) std::cout<<"mark elements"<<std::endl;
+ Dune::PDELab::mark_grid(grid,z0,eta_refine,eta_coarsen, Rminlevel,Rmaxlevel,verbosity);
+ }
+};
diff --git a/dune/phasefield/debug/vector_debug.hh b/dune/phasefield/debug/vector_debug.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1a81695d7027bac516aea66d122f2506ca1b3a74
--- /dev/null
+++ b/dune/phasefield/debug/vector_debug.hh
@@ -0,0 +1,36 @@
+#pragma once
+#include "../myincludes.hh"
+
+#include <string>
+
+
+
+
+
+
+
+template<typename V>
+int calculateVectorSize(const V& v)
+{
+ /*int i=0;
+ for(auto it=v.begin(); it!=v.end(); ++it)
+ {
+ ++i;
+ }
+ return i;*/
+ return v.flatsize();
+}
+
+template<typename V>
+void printVector(const V& v, const int maxnum, std::string Name = "")
+{
+ int size = calculateVectorSize(v);
+ int maxiter = std::min(maxnum, size);
+ std::cout << "Printing vector '" << Name << "' (" << maxnum << " of " << size << "):" << std::endl;
+ auto it = v.begin();
+ for(int i=0; i<maxiter; ++i)
+ {
+ std::cout<< *it <<std::endl;
+ ++it;
+ }
+}
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
new file mode 100644
index 0000000000000000000000000000000000000000..52fb625efb8e1416ee6558a3f6d613710c853d94
--- /dev/null
+++ b/dune/phasefield/initial_fn.hh
@@ -0,0 +1,58 @@
+
+#pragma once
+#include "myincludes.hh"
+#include <cmath>
+
+
+template<typename GV, typename RF, std::size_t dim_range, typename Params>
+class ZeroInitialVector :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range>,ZeroInitialVector<GV,RF,dim_range, Params> >
+{
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range> Traits;
+
+ ZeroInitialVector(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits, ZeroInitialVector>(gv) {}
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ y=0;
+ }
+};
+
+//______________________________________________________________________________________________________________________
+
+template<typename GV, typename RF, typename Params>
+class ZeroInitial
+ : public ZeroInitialVector<GV,RF,1,Params>
+{
+public:
+ ZeroInitial(const GV & gv, const Params & params) : ZeroInitialVector<GV,RF,1,Params>(gv,params) {}
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename GV, typename RF, typename Params>
+class PhiInitial_RayleighTaylor :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_RayleighTaylor<GV,RF,Params> >
+{
+ RF eps;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_RayleighTaylor(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_RayleighTaylor<GV,RF,Params> > (gv)
+ {
+ eps = params.eps;
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ y = 0.5*(1-tanh(1/eps/sqrt(2)*(x[1]-0.75-0.01*cos(x[0]*2*3.1415))));
+ }
+};
+
+//_________________________________________________________________________________________________________________________
diff --git a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
index 3d7e907dd8a9042f7238dbd20e8378a7c4311d47..bd1ac1f714428794212fc43281ded731c1300e37 100644
--- a/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
+++ b/dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh
@@ -1,6 +1,6 @@
#pragma once
-#include "myincludes.hh"
+#include "../myincludes.hh"
#include <dune/istl/io.hh>
#include <dune/istl/operators.hh>
diff --git a/src/test_2f0s/acfs_fem_ac_estimator.hh b/src/test_2f0s/acfs_fem_ac_estimator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..af958ce5ee17fe30d0ab5d25b662c91650346ab3
--- /dev/null
+++ b/src/test_2f0s/acfs_fem_ac_estimator.hh
@@ -0,0 +1,64 @@
+#pragma once
+
+#include "myincludes.hh"
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+
+
+template<typename CHGFS>
+class ACFS_FEM_AC_Estimator :
+ public Dune::PDELab::LocalOperatorDefaultFlags
+{
+private:
+ typedef typename CHGFS::template Child<0>::Type FEM;
+ typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
+ Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+ enum {dim=LocalBasis::Traits::dimDomain};
+
+public:
+ enum {doPatternVolume=false};
+ //enum {doLambdaVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ ACFS_FEM_AC_Estimator()
+ {
+ }
+
+
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+ {
+ using namespace Dune::Indices;
+ const auto& pfspace = child(lfsu,_0);
+ typedef decltype(makeZeroBasisFieldValue(pfspace)) RF;
+
+ auto geo = eg.geometry();
+ const auto& rElement = Dune::ReferenceElements<double, dim>::general(geo.type());
+ RF umax=0.0;
+ RF umin=0.0;
+ // Dune::FieldVector<RF, dim> xg = eg.geometry().global( rElement.position(0, dim) );
+ for (int i=0; i<geo.corners(); i++)
+ {
+ RF u=0.0;
+ auto& phihat = cache.evaluateFunction(rElement.position(i, dim),pfspace.finiteElement().localBasis());
+ for (size_t j=0; j<pfspace.size(); j++)
+ u += x(pfspace,j)*phihat[j];
+
+ if(i==0)
+ {
+ umax = u;
+ umin = u;
+ }
+ else
+ {
+ if(u>umax) umax = u;
+ if(u<umin) umin = u;
+ }
+ //xg = eg.geometry().global( rElement.position(i, dim) );
+ }
+ r.accumulate(lfsv,0,(umax-umin));//(xg[0]>0.25&&xg[0]<0.75)? (umax-umin):0 );
+ }
+
+};
diff --git a/src/test_2f0s/chns.ini b/src/test_2f0s/chns.ini
index 22ba3e00c948c03a87b6617dd576afc81953dc8a..504f7802e05693e5a0fb97f83057f693afe1fec1 100644
--- a/src/test_2f0s/chns.ini
+++ b/src/test_2f0s/chns.ini
@@ -24,12 +24,14 @@ filename = grids/rectangle_3k_per.msh
#filename = grids/square_2k_del_per.msh
-[temporal]
-time = 30
-taumax = 0.01
-tau = 0.015
-taumin = 0.0001
-saveintervall = 0.05
+[Time]
+tMax = 30
+dt = 0.015
+#dtMax = 0.015
+#dtMin = 0.015
+#saveintervall = 0.5
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
[problem]
eps=0.03
@@ -42,3 +44,9 @@ delta=0.03
u=1
vmax=0
vmaxtime=5
+
+[Refinement]
+etaRefine = 0.15
+etaCoarsen = 0.05
+maxRefineLevel = 2
+minRefineLevel = 0
diff --git a/src/test_2f0s/chns_driver.hh b/src/test_2f0s/chns_driver.hh
index 29cddcddd98eb3aa376aec0e76e8f8cc60c35b39..e193058a578e5b2b4bd0252aa15b1ea877ff7f07 100644
--- a/src/test_2f0s/chns_driver.hh
+++ b/src/test_2f0s/chns_driver.hh
@@ -11,7 +11,7 @@
//#include "acfs_fem_allencahn.hh"
//#include "acfs_fem_cahnhilliard.hh"
//#include "acfs_fem_navierstokes.hh"
-//#include "acfs_fem_ac_estimator.hh"
+#include "acfs_fem_ac_estimator.hh"
//#include "threephase_fem_cahnhilliard.hh"
//#include "threephase_fem_navierstokes.hh"
//#include "threephase_dgf_coloring.hh"
@@ -22,7 +22,9 @@
#include <dune/phasefield/acfs_ibvp.hh>
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
-
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
//template <class A>
//using MySolver = Dune::RestartedGMResSolver<A, A, A>;
@@ -41,191 +43,76 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
const int dim = GV::dimension;
typedef double RF; // type for computations
- RF T = ptree.get("temporal.time",(RF)1.0);
RF initialvmax = ptree.get("initial.vmax",(RF)1);
RF vmaxtime = ptree.get("initial.vmaxtime",(RF)1);
- RF time = 0.0;
- RF lastdatasave = 0.0;
- RF datasaveintervall = ptree.get("temporal.saveintervall",(RF)1.0);
- RF dt_min = 1e-8;
- bool savealways = ptree.get("output.savealways",(RF)1);
-
- Problem<RF> problem(ptree.get("temporal.tau",(RF)0.1),ptree);
- problem.setTime(time);
-
-//_________________________________________________________________________________________________________________________
-//Initial Values
- //typedef ZeroScalarFunction<GV,RF> InitialPressure;
- //InitialPressure initialPressure(gv);
- typedef TUPressure<GV,RF,1> InitialPressure;
- InitialPressure initialPressure(gv);
-
- typedef IniVelocity<GV,RF,2> InitialVelocity;
- InitialVelocity initialVelocity(gv,initialvmax,vmaxtime);
- typedef Dune::PDELab::CompositeGridFunction<InitialVelocity,InitialPressure> InitialNS;
- InitialNS initialNS(initialVelocity,initialPressure);
+ TimestepManager timestepManager(ptree.sub("Time"),2);
- typedef PhiZero<GV,RF,1> InitialPhi;
- InitialPhi initialPhi(gv,0.5,0.15,problem.eps);
- typedef ThreephasePhiZero<GV,RF,1> InitialThreephasePhi;
- InitialThreephasePhi initialPhiOne(gv,0, 0.5,0.15,problem.eps);
- InitialThreephasePhi initialPhiTwo(gv,1, 0.5,0.15,problem.eps);
- typedef ZeroScalarFunction<GV,RF> InitialMu;
- InitialMu initialMu(gv);
-
- typedef Dune::PDELab::CompositeGridFunction<InitialPhi,InitialMu> InitialCH;
- InitialCH initialCH(initialPhi,initialMu);
+ Problem<RF> problem(ptree.get("temporal.tau",(RF)0.1),ptree);
+ problem.setTime(timestepManager.t);
//_________________________________________________________________________________________________________________________
// Make grid function spaces
- typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
- //typedef Dune::PDELab::NoConstraints Unconstraint;
- typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend; //Dune::PDELab::ISTL::Blocking::none
- typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
- //typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::bcrs> VB_Block2;
- using OrderingTagL = Dune::PDELab::LexicographicOrderingTag;
- using OrderingTagB = Dune::PDELab::EntityBlockedOrderingTag;
-
- ConstraintsAssembler conp, conphi, conmu;
-
- typedef Dune::PDELab::VectorGridFunctionSpace
- <ES,V_FEM,dim,VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
- V_GFS vGfs(es,vFem);
- vGfs.name("velocity");
-
- /*typedef Dune::PDELab::GridFunctionSpace<ES, V_FEM, ConstraintsAssembler, VectorBackend> VK_GFS;
- VK_GFS v1Gfs(es, vFem, conv1);
- VK_GFS v2Gfs(es, vFem, conv2);
-
- typedef Dune::PDELab::PowerGridFunctionSpace<VK_GFS, dim, VectorBackend, OrderingTagL> V_GFS;
- V_GFS vGfs(v1Gfs, v2Gfs);
- vGfs.name("velocity");*/
-
- typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
- P_GFS pGfs(es,pFem, conp);
- pGfs.name("pressure");
-
- typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> PHI_GFS;
- PHI_GFS phiGfs(es,phiFem, conphi);
- phiGfs.name("phasefield");
-
- typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> MU_GFS;
- MU_GFS muGfs(es,phiFem, conmu);
- muGfs.name("potential");
-
- typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS_GFS;
- NS_GFS nsGfs(vGfs, pGfs);
- typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS> CH_GFS;
- CH_GFS chGfs(phiGfs, muGfs);
-
-
- //typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, NS_GFS, CH_GFS> GFS;
- //GFS gfs(nsGfs, chGfs);
+ typedef GFS_2f0s<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+ GFS gfs(es,vFem,pFem,phiFem);
+ typedef typename GFS::NS NS_GFS;
+ typedef typename GFS::CH CH_GFS;
//_________________________________________________________________________________________________________________________
//Constraints
- //typedef Dune::PDELab::EmptyTransformation CC;
- typedef typename NS_GFS::template ConstraintsContainer<RF>::Type NSC;
- NSC nscg;
- nscg.clear();
-
- typedef typename CH_GFS::template ConstraintsContainer<RF>::Type CHC;
- CHC chcg;
-
- // create Taylor-Hood constraints (no-slip, no penetration)
- typedef Dune::PDELab::DirichletConstraintsParameters AlwaysDirichlet;
- typedef Dune::PDELab::NoDirichletConstraintsParameters NeverDirichlet;
- typedef Dune::PDELab::PowerConstraintsParameters<AlwaysDirichlet,dim> VelocityConstraints;
- typedef Dune::PDELab::CompositeConstraintsParameters<VelocityConstraints, SinglePointPressureBC> NSConstraints;
-
- AlwaysDirichlet scalarvelocity_constraints;
- SinglePointPressureBC pressure_constraints;
- VelocityConstraints velocity_constraints(scalarvelocity_constraints);
- NSConstraints nsConstraints(velocity_constraints,pressure_constraints);
-
- NeverDirichlet pf_constraints;
- NeverDirichlet mu_constraints;
-
- typedef Dune::PDELab::CompositeConstraintsParameters<NeverDirichlet, NeverDirichlet> CHChonstraints;
- CHChonstraints chConstraints(pf_constraints,mu_constraints);
-
-
-
- Dune::PDELab::constraints(nsConstraints,nsGfs,nscg, 0);
- Dune::PDELab::constraints(chConstraints,chGfs,chcg, 0);
+ typedef Bdry_2f0s<dim, Dune::PDELab::DirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
+ Bdry bdry;
+ gfs.updateConstraintsContainer(bdry);
//_________________________________________________________________________________________________________________________
//Coefficient vectors
using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
- CH_V chV(chGfs);
- NS_V nsV(nsGfs);
- chV=1.1245;
- nsV=2;
+ CH_V chV(gfs.ch);
+ NS_V nsV(gfs.ns);
CH_V chVOld = chV;
NS_V nsVOld = nsV;
- CH_V chInterpolant = chV;
- NS_V nsInterpolant = nsV;
-//_________________________________________________________________________________________________________________________
+ //_________________________________________________________________________________________________________________________
+ //Initial Values
+
+ typedef Ini_2f0s< ZeroInitialVector<GV,RF,dim,Problem<RF>>,
+ ZeroInitial<GV,RF,Problem<RF>>,
+ PhiInitial_RayleighTaylor<GV,RF,Problem<RF>>,
+ ZeroInitial<GV,RF,Problem<RF>> > Initial;
+ Initial initial(gv,problem);
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+ //_________________________________________________________________________________________________________________________
// LocalOperators
typedef CHNS_FEM_NavierStokes<Problem<RF>, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
- NS_LOP nsLop(problem, chGfs, chV, chVOld, nsGfs, nsVOld);
+ NS_LOP nsLop(problem, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
typedef CHNS_FEM_CahnHilliard<Problem<RF>, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
- CH_LOP chLop(problem, chGfs, chVOld, nsGfs, nsV);
+ CH_LOP chLop(problem, gfs.ch, chVOld, gfs.ns, nsV);
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
-
-
-
- //_________________________________________________________________________________________________________________________
- // Interpolate Initial Data
- //Dune::PDELab::interpolate(initialFn,nsGfs,nsVOld);
-
- Dune::PDELab::interpolate(initialCH,chGfs,chV);
- Dune::PDELab::interpolate(initialNS,nsGfs,nsV);
-
-
-
//_________________________________________________________________________________________________________________________
// Construct discrete Grid Functions
- typedef typename Dune::PDELab::GridFunctionSubSpace<NS_GFS,Dune::TypeTree::TreePath<0> > VelocitySubGFS;
- VelocitySubGFS velocitySubGfs(nsGfs);
- typedef typename Dune::PDELab::GridFunctionSubSpace<NS_GFS,Dune::TypeTree::TreePath<1> > PressureSubGFS;
- PressureSubGFS pressureSubGfs(nsGfs);
- typedef Dune::PDELab::VectorDiscreteGridFunction<VelocitySubGFS,NS_V> VDGF;
- VDGF vDgf(velocitySubGfs,nsV);
- typedef Dune::PDELab::DiscreteGridFunction<PressureSubGFS,NS_V> PDGF;
- PDGF pDgf(pressureSubGfs,nsV);
-
-
- typedef typename Dune::PDELab::GridFunctionSubSpace<CH_GFS,Dune::TypeTree::TreePath<0> > PhasefieldSubGFS;
- PhasefieldSubGFS phasefieldSubGfs(chGfs);
-
- typedef typename Dune::PDELab::GridFunctionSubSpace<CH_GFS,Dune::TypeTree::TreePath<1> > PotentialSubGFS;
- PotentialSubGFS potentialSubGfs(chGfs);
-
- typedef Dune::PDELab::DiscreteGridFunction<PotentialSubGFS,CH_V> MuDGF;
- MuDGF muDgf(potentialSubGfs,chV);
-
-
- typedef Dune::PDELab::DiscreteGridFunction<PhasefieldSubGFS,CH_V> PhiDGF;
- PhiDGF phiDgf(phasefieldSubGfs,chV);
+ typedef DGF_2f0s<GFS, NS_V, CH_V> DGF;
+ DGF dgf(gfs, nsV, chV);
//_________________________________________________________________________________________________________________________
// prepare VTK writer
@@ -241,71 +128,38 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
}
auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vDgf,"v"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
+ dgf.addToVTKWriter(vtkwriter);
+ vtkwriter.write(timestepManager.t,Dune::VTK::ascii);
//_________________________________________________________________________________________________________________________
- int i=0;
- for(auto it = chV.begin(); it!=chV.end(); ++it)
- {
- ++i;
- if(i<100) std::cout<< *it <<std::endl;
- }
- std::cout << i << std::endl;
- i=0;
- for(auto it = nsV.begin(); it!=nsV.end(); ++it)
- {
- ++i;
- if(i<100) std::cout<< *it <<std::endl;
- }
- std::cout << i << std::endl;
-
- Dune::PDELab::set_constrained_dofs(nscg,0.0,nsV);
- //Dune::PDELab::set_constrained_dofs(chcg,0.0,chV);
-
- vtkwriter.write(time,Dune::VTK::ascii);
+ printVector(chV,50,"chV");
+ printVector(nsV,50,"nsV");
- //_________________________________________________________________________________________________________________________
+ //____________________________________________________________________________________________________________________
//int i_it=0;
chVOld = chV;
nsVOld = nsV;
- while (time < T - dt_min*0.5)
+ while (!timestepManager.isFinished())
{
- problem.setTime(time+problem.dt);
+ problem.setTime(timestepManager.t);
// Grid Operators
- typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,NSC,NSC> NS_GO;
- NS_GO nsGo(nsGfs,nscg,nsGfs,nscg,nsLop,mbe);
+ typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
+ NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
- typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,CHC,CHC> CH_GO;
- CH_GO chGo(chGfs,chcg,chGfs,chcg,chLop,mbe);
+ typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
+ CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
// Linear solver
- //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 CH_LS; CH_LS ls_ch(200,5000,1);
- //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_SSORk<CH_GO> CH_LS; CH_LS ls_ch(chGo,5000,3,1);
-
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
<CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,3,1);
-
- //typedef Dune::PDELab::ISTLBackend_SEQ_GMRES_ILU0 NS_LS; NS_LS ls_ns(200,5000,1);
- //typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_ILU0 NS_LS; NS_LS ls_ns(5000,1);
-
- //AT LEAST WORKS:
- //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_NOPREC<NS_GFS> NS_LS; NS_LS ls_ns(nsGfs,5000,1);
-
- //typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_Jacobi<NS_GFS> NS_LS; NS_LS ls_ns(nsGfs,5000,1);
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
<NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
- //typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_Jac LS;
- //typedef Dune::PDELab::ISTLBackend_SEQ_SuperLU LS;
- //typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LS;
- // Solve (possibly) nonlinear problem
+
typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
PDESolver1 ns_newton(nsGo,nsV,ls_ns);
ns_newton.setReassembleThreshold(0.0);
@@ -327,118 +181,51 @@ void chns_driver (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem,
//_________________________________________________________________________________________________________________________
// do time step
- initialNS.setTime(time);
- Dune::PDELab::interpolate(initialNS,nsGfs,nsInterpolant);
- Dune::PDELab::interpolate(initialCH,chGfs,chInterpolant);
- Dune::PDELab::copy_constrained_dofs(nscg, nsInterpolant, nsV);
- Dune::PDELab::copy_constrained_dofs(chcg, chInterpolant, chV );
- bool EverythingOk = true;
- try
+
+
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
+
+ int result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+ if(result <0)
{
- if (helper.rank()==0)
- std::cout << "=== Begin newtonapply CH:" << std::endl;
- ch_newton.apply();
- //vtkwriter.write(time+5e-6,Dune::VTK::ascii);
- if (helper.rank()==0)
- std::cout << "=== Begin newtonapply NS:" << std::endl;
- ns_newton.apply();
+ return;
}
- catch(...)
+ else if(result == 0)
{
- //vtkwriter.write(time+6e-6,Dune::VTK::ascii);
- if (helper.rank()==0)
- std::cout << "Reducing Timestep to "<< problem.dt/2 << std::endl;
- problem.dt = problem.dt/2;
- EverythingOk = false;
- //time+=1e-5;
- //NS_V nsRes(nsGfs);
- //nsGo.residual(nsV,nsRes);
- //nsV = nsRes;
- //vtkwriter.write(time,Dune::VTK::ascii);
nsV = nsVOld;
chV = chVOld;
- if(problem.dt<ptree.get("temporal.taumin",(RF)0.0001))
- {
- if (helper.rank()==0)
- std::cout << "Timestep too small, aborting :(" << std::endl;
- return;
- }
+ continue;
}
- if(EverythingOk)
- {
- time += problem.dt;
- problem.dt = fmin(1.03*problem.dt, ptree.get("temporal.taumax",(RF)1000) );
- if (helper.rank()==0)
- std::cout << "New Timestep "<< problem.dt << std::endl;
- if( (time>lastdatasave+datasaveintervall) || savealways)
- {
- lastdatasave = time;
- vtkwriter.write(time,Dune::VTK::ascii);
- }
- chVOld = chV;
- nsVOld = nsV;
- }
+ chVOld = chV;
+ nsVOld = nsV;
//_________________________________________________________________________________________________________________________
-/*
- //Estimator GFS
- typedef Dune::PDELab::P0LocalFiniteElementMap<RF,RF,dim> P0FEM;
- P0FEM p0fem(Dune::GeometryTypes::simplex(dim));
- typedef Dune::PDELab::NoConstraints NCON;
- typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM,NCON,VectorBackend> P0GFS;
- P0GFS p0Gfs(gv,p0fem);
- //Estimator vector
- using Z0 = Dune::PDELab::Backend::Vector<P0GFS,RF>;
- Z0 z0(p0Gfs,0.0);
- typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
- Z0DGF z0dgf(p0Gfs,z0);
- //Estimator Operator
- typedef ACFS_FEM_AC_Estimator<Problem<RF>,CH_GFS> EST_LOP;
- EST_LOP estlop(problem);
- typedef Dune::PDELab::EmptyTransformation NCC;
- typedef Dune::PDELab::GridOperator<CH_GFS,P0GFS,EST_LOP,MBE,RF,RF,RF,NCC,NCC> ESTGO;
- ESTGO estgo(chGfs,p0Gfs,estlop,mbe);
- estgo.residual(chV,z0);
+
// mark elements for refinement
- std::cout << "mark elements" << std::endl;
- RF eta_refine = ptree.get("fem.etarefine",(double)0.15);
- RF eta_coarsen = ptree.get("fem.etarecoarsen",(double)0.0);
- int Rmaxlevel = ptree.get("fem.maxlevel",(double)3);
- int Rminlevel = ptree.get("fem.minlevel",(double)2);
- Dune::PDELab::mark_grid(grid,z0,eta_refine,eta_coarsen, Rminlevel,Rmaxlevel,2);
+ typedef RefinementIndicator<GV,RF,ACFS_FEM_AC_Estimator<CH_GFS>,CH_GFS> RIND;
+ RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+ rInd.calculateIndicator(chV);
+ rInd.markGrid(grid,2);
- //znew = z;
- //znew = 1.0;
- std::cout << "Size Before:" << chGfs.globalSize() << std::endl;
+ printVector(chV,3,"chV");
+ printVector(nsV,3,"nsV");
// do refinement
std::cout << "adapt grid and solution" << std::endl;
- auto chTransfer = transferSolutions(chGfs,2,chV,chVOld);
- auto nsTransfer = transferSolutions(nsGfs,2,nsV,nsVOld);
+ auto chTransfer = transferSolutions(gfs.ch,2,chV,chVOld);
+ auto nsTransfer = transferSolutions(gfs.ns,2,nsV,nsVOld);
Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
- //Dune::PDELab::adapt_grid(grid, chGfs, chV, 2*2);
std::cout << "adapt grid and solution ended" << std::endl;
- //znew = z;
- //znew = 1.0;
- std::cout << "Size After:" << chGfs.globalSize() << std::endl;
+
+ printVector(chV,3,"chV");
+ printVector(nsV,3,"nsV");
+
// recompute constraints
- Dune::PDELab::constraints(nsConstraints,nsGfs,nscg, 0);
- Dune::PDELab::constraints(chConstraints,chGfs,chcg, 0);
- Dune::PDELab::interpolate(initialNS,nsGfs,nsVOld);
- Dune::PDELab::interpolate(initialCH,chGfs,chVOld);
- if(i_it>5)
- {
- Dune::PDELab::copy_constrained_dofs(nscg, nsV , nsVOld );
- Dune::PDELab::copy_constrained_dofs(chcg, chV , chVOld );
- }
- else
- {
- chV=chVOld;
- nsV=nsVOld;
- }
- i_it++;
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+
- vtkwriter.write(time+1e-6,Dune::VTK::ascii);*/
+ vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
}
}