diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index 5fe527fe157cba2f6be62db44ba60e3bc1f30151..41c6ad9f81d4b760af464ecf4cc387aca9908975 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -463,7 +463,7 @@ void my_adapt_grid (Grid& grid, GV& gv, GFS1& gfs1, X1& x1, X1& x1b, GFS2& gfs2,
}
-
+#define BLOCKED
//_________________________________________________________________________________________________________________________
template<typename GV, typename RF, typename EST_LOP, typename EST_GFS>
@@ -474,14 +474,23 @@ public:
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;
+
+#ifndef BLOCKED
+ typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM,NCON,Dune::PDELab::ISTL::VectorBackend<>> P0GFS;
+#else
+ using BlockedVBE = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::bcrs>;
+ typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
+ typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM,NCON,Dune::PDELab::ISTL::VectorBackend<>> P1GFS;
+ typedef Dune::PDELab::CompositeGridFunctionSpace<BlockedVBE, OrderingTagL, P1GFS> P0GFS;
+ P1GFS p1Gfs;
+#endif
P0GFS p0Gfs;
//Estimator vector
using Z0 = Dune::PDELab::Backend::Vector<P0GFS,RF>;
Z0 z0;
- typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
- Z0DGF z0dgf;
+// typedef Dune::PDELab::DiscreteGridFunction<P0GFS,Z0> Z0DGF;
+// Z0DGF z0dgf;
//Estimator Operator
EST_LOP estlop;
@@ -497,9 +506,14 @@ public:
RefinementIndicator(const GV& gv, const EST_GFS chGfs, const Dune::ParameterTree & param) :
p0fem(Dune::GeometryTypes::simplex(GV::dimension)),
+#ifndef BLOCKED
p0Gfs(gv,p0fem),
+#else
+ p1Gfs(gv,p0fem),
+ p0Gfs(p1Gfs),
+#endif
z0(p0Gfs,0.0),
- z0dgf(p0Gfs,z0),
+// z0dgf(p0Gfs,z0),
estlop(),
estgo(chGfs,p0Gfs,estlop,MBE(1))
{
@@ -516,7 +530,7 @@ public:
if(writeOutput)
{
auto vtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,Dune::refinementLevels(0));
- vtkwriter->addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Z0DGF> >(z0dgf,"indicator"));
+// vtkwriter->addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Z0DGF> >(z0dgf,"indicator"));
vtkwriter->write("RefinementIndicator",Dune::VTK::ascii);
}
return z0;
diff --git a/dune/phasefield/ff_chns.hh b/dune/phasefield/ff_chns.hh
index 4e1d601b5fbec9209d69b07066adb5671245fbf8..e1c5442e753e7e6aa95fa412e66001dd0ea85768 100644
--- a/dune/phasefield/ff_chns.hh
+++ b/dune/phasefield/ff_chns.hh
@@ -18,13 +18,14 @@ public:
typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
+ typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::bcrs> VB_BCRS;
typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
ConstraintsAssembler conp, conphi, conmu;
- typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
- VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+// typedef Dune::PDELab::VectorGridFunctionSpace<ES, V_FEM, ES::GridView::dimension, VectorBackend, VectorBackend, ConstraintsAssembler> V_GFS;
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES, V_FEM, ES::GridView::dimension, VectorBackend, VectorBackend, ConstraintsAssembler, OrderingTagB> V_GFS;
V_GFS vGfs;
typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
@@ -39,7 +40,9 @@ public:
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;
+// typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS> CH;
+// typedef Dune::PDELab::CompositeGridFunctionSpace<VB_BCRS, OrderingTagB, PHI_GFS, MU_GFS> CH;
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagB, PHI_GFS, MU_GFS> CH;
CH ch;
typedef typename NS::template ConstraintsContainer<RF>::Type NS_CC;
@@ -68,8 +71,72 @@ public:
template<typename Bdry>
void updateConstraintsContainer(const Bdry& bdry)
{
- Dune::PDELab::constraints(bdry.bdryNS,ns,nsCC, 0);
- Dune::PDELab::constraints(bdry.bdryCH,ch,chCC, 0);
+ Dune::PDELab::constraints(bdry.bdryNS, ns, nsCC, 0);
+ Dune::PDELab::constraints(bdry.bdryCH, ch, chCC, 0);
+ }
+};
+
+template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+class GFS_ff_chns_blocked
+{
+public:
+ typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
+ typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
+
+ // define function space options
+ using FlatVBE = Dune::PDELab::ISTL::VectorBackend<>;
+ using BlockedVBE = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::bcrs>;
+
+ ConstraintsAssembler conp, conphi, conmu;
+
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES, V_FEM, ES::GridView::dimension, FlatVBE, FlatVBE, ConstraintsAssembler, OrderingTagB> V_GFS;
+ V_GFS vGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, FlatVBE> P_GFS;
+ P_GFS pGfs;
+
+ // make function space tree
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, FlatVBE> PHI_GFS;
+ PHI_GFS phiGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, FlatVBE> MU_GFS;
+ MU_GFS muGfs;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<BlockedVBE, OrderingTagL, V_GFS, P_GFS> NS;
+ NS ns;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<BlockedVBE, OrderingTagL, PHI_GFS, MU_GFS> CH;
+ CH ch;
+
+ typedef typename NS::template ConstraintsContainer<RF>::Type NS_CC;
+ NS_CC nsCC;
+
+ typedef typename CH::template ConstraintsContainer<RF>::Type CH_CC;
+ CH_CC chCC;
+
+ GFS_ff_chns_blocked(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.name("phasefield");
+ 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);
}
};
@@ -79,14 +146,14 @@ class GFS_ff_chns_monolithic
public:
typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
- typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
+ typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::bcrs> VB_Block;
typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
ConstraintsAssembler conp, conphi, conmu;
- typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
- VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+// typedef Dune::PDELab::VectorGridFunctionSpace<ES, V_FEM, ES::GridView::dimension, VectorBackend, VectorBackend, ConstraintsAssembler> V_GFS;
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES, V_FEM, ES::GridView::dimension, VectorBackend, VectorBackend, ConstraintsAssembler, OrderingTagB> V_GFS;
V_GFS vGfs;
typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
@@ -134,6 +201,74 @@ public:
}
};
+template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+class GFS_ff_chns_monolithic_blocked
+{
+public:
+ typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
+ typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
+ typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
+ typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
+
+ ConstraintsAssembler conp, conphi, conmu;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, V_FEM, ConstraintsAssembler, VectorBackend> Vx_GFS;
+ Vx_GFS vxGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, V_FEM, ConstraintsAssembler, VectorBackend> Vy_GFS;
+ Vy_GFS vyGfs;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagB, Vx_GFS, Vy_GFS> V_GFS;
+ V_GFS vGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> PHI_GFS;
+ PHI_GFS phiGfs;
+
+ 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;
+ NS ns;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagB, PHI_GFS, MU_GFS> CH;
+ CH ch;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, CH, NS> CHNS;
+ CHNS chns;
+
+ typedef typename CHNS::template ConstraintsContainer<RF>::Type CHNS_CC;
+ CHNS_CC chnsCC;
+
+ GFS_ff_chns_monolithic_blocked(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem) :
+ conp(), conphi(), conmu(),
+ vxGfs(es, vFem),
+ vyGfs(es, vFem),
+ vGfs(vxGfs, vyGfs),
+ pGfs(es, pFem, conp),
+ phiGfs(es, phiFem, conphi),
+ muGfs(es, phiFem, conmu),
+ ns(vGfs, pGfs),
+ ch(phiGfs, muGfs),
+ chns(ch, ns),
+ chnsCC()
+ {
+ vGfs.name("velocity");
+ pGfs.name("pressure");
+ phiGfs.name("phasefield");
+ muGfs.name("potential");
+ }
+
+ template<typename Bdry>
+ void updateConstraintsContainer(const Bdry& bdry)
+ {
+ Dune::PDELab::constraints(bdry.bdryCHNS, chns, chnsCC, 0);
+ }
+};
+
//_________________________________________________________________________________________________________________________
@@ -298,7 +433,7 @@ public:
BdryP bdryP;
BdryPhi bdryPhi;
BdryMu bdryMu;
- typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> BdryV;
+ typedef Dune::PDELab::PowerConstraintsParameters<BdryV_Scalar,dim> BdryV;
typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
BdryV bdryV;
@@ -357,7 +492,7 @@ public:
BdryP bdryP;
BdryPhi bdryPhi;
BdryMu bdryMu;
- typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2> BdryV;
typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
BdryV bdryV;
diff --git a/src/example_2p_cahnhilliard/2p_driver.hh b/src/example_2p_cahnhilliard/2p_driver.hh
index bec949ed926b1ba65af6a325a5b946aca837dff7..acd01aa9a0a9d1561d7d764cc3ac3075424573c0 100644
--- a/src/example_2p_cahnhilliard/2p_driver.hh
+++ b/src/example_2p_cahnhilliard/2p_driver.hh
@@ -30,7 +30,7 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::Paramete
typedef double RF;
// We use the Paramter File for the Cahn--Hilliard Navier--Stokes Model,
- //as it includes all relevant parameters of the Cahn--Hilliard Model
+ // as it includes all relevant parameters of the Cahn--Hilliard Model
typedef Params_ff< RF, DoubleWell_poly<RF> , Mobility_Const<RF> > Parameters;
Parameters param(ptree);
@@ -41,14 +41,15 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::Paramete
// Again, we use the Cahn--Hilliard Navier--Stokes Model, in operator splitting form the Cahn--Hilliard grid function space
// is exactly the space we need.
- typedef GFS_ff_chns<RF, ES, PHI_FEM, PHI_FEM, PHI_FEM> GFS;
+// typedef GFS_ff_chns<RF, ES, PHI_FEM, PHI_FEM, PHI_FEM> GFS;
+ typedef GFS_ff_chns_blocked<RF, ES, PHI_FEM, PHI_FEM, PHI_FEM> GFS;
GFS gfs(es,phiFem,phiFem,phiFem);
typedef typename GFS::NS NS_GFS;
typedef typename GFS::CH CH_GFS;
//_________________________________________________________________________________________________________________________
- //Constraints
+ // Constraints
typedef Bdry_ff_chns<2, Dune::PDELab::DirichletConstraintsParameters, // v_x, v_y
Dune::PDELab::NoDirichletConstraintsParameters, // p
@@ -59,7 +60,7 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::Paramete
gfs.updateConstraintsContainer(bdry);
//_________________________________________________________________________________________________________________________
- //Coefficient vectors
+ // Coefficient vectors
using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
@@ -67,11 +68,11 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::Paramete
NS_V nsV(gfs.ns);
//_________________________________________________________________________________________________________________________
- //Initial Values
+ // Initial Values
typedef Ini_ff_chns< ZeroInitial<GV,RF,Parameters>, //v
ZeroInitial<GV,RF,Parameters>, //p
- PhiInitial_SpinodalDecomposition<GV,RF,Parameters>, //phi
+ PhiInitial_Sphere<GV,RF,Parameters>, //phi
ZeroInitial<GV,RF,Parameters>> Initial; //mu
Initial initial(gv,param);
Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
@@ -87,6 +88,8 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::Paramete
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+// typedef Dune::PDELab::Backend::Matrix<Dune::PDELab::ISTL::BCRSMatrixBackend<>,CH_V,CH_V,double> MBE;
+// MBE mbe();
//_________________________________________________________________________________________________________________________
@@ -144,9 +147,11 @@ void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::Paramete
CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
// Linear solver
- typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
- <CH_GO, Dune::SeqILU, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,600,200,3,1);
-
+// typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+// <CH_GO, Dune::SeqILU, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,600,200,3,1);
+ typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_Richardson CH_LS;
+ CH_LS ls_ch(10000, 1);
+
typedef Dune::PDELab::NewtonMethod<CH_GO, CH_LS> PDESolver2;
PDESolver2 ch_newton(chGo, ls_ch);
ch_newton.setParameters(ptree.sub("Phasefield_Newton"));