diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index f1b033cd1ecd93b5db8d60017439d59bc6b016f2..dd4b7e1340857f788438626cc16ebe3cda01b3ce 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -50,8 +50,7 @@ public:
Number eval( const Number& pf1, const Number& pf2, const Number& pf3) const
{
- std::cout << "Mobility_Degen for three phases not properly implemented yet" << std::endl;
- return mobility;
+ return std::max(mobility*(1-pf3)*(1-pf3) + mobility*pf1*pf2*pf3*9,mobility_reg);
}
};
@@ -362,6 +361,19 @@ public:
{
}
};
+
+template <typename Number, typename TW, typename Mobility, typename Reaction, typename VelDissipation>
+class Params_summerschool_3p : public Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>
+{
+public:
+ Number CellHeight;
+
+ Params_summerschool_3p(const Dune::ParameterTree& ptree) :
+ Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>(ptree),
+ CellHeight(ptree.get("domain.CellHeight",0.1))
+ {
+ }
+};
//
// template<typename Number>
// class Problem
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index a1f760fb21e8e4b002c4045c9afc24901a6f560e..605fde567328393a1a3d58545df6da86083fe86e 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -387,10 +387,10 @@ public:
}
};
-template<typename GV, typename RF, typename Params>
+template<typename GV, typename RF, typename Params, const int phase = -1>
class PhiInitial_SummerSchoolPaper :
public Dune::PDELab::AnalyticGridFunctionBase<
- Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_SummerSchoolPaper<GV,RF,Params> >
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_SummerSchoolPaper<GV,RF,Params,phase> >
{
public:
RF eps;
@@ -400,7 +400,8 @@ public:
private:
inline static RF PxToCoord_x(const RF px_x)
{
- return 0.5+(px_x - 150)/((RF)299);
+ //return 0.5+(px_x - 150)/((RF)299);
+ return 0.5+(px_x - 150)/((RF)299)-3.5;
}
inline static RF PxToCoord_y(const RF px_y)
{
@@ -437,37 +438,42 @@ public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
PhiInitial_SummerSchoolPaper(const GV & gv, const Params & params) :
- Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_SummerSchoolPaper<GV,RF,Params> > (gv)
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_SummerSchoolPaper<GV,RF,Params,phase> > (gv)
{
eps = params.eps;
radius = params.initial.ptree.get("Radius",0.2);
-/*
+
//Experiment 1
- xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75));
+ /* xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75));
xValues.push_back(PxToCoord_x(317)); yValues.push_back(PxToCoord_y(85));
xValues.push_back(PxToCoord_x(364)); yValues.push_back(PxToCoord_y(110));
xValues.push_back(PxToCoord_x(556)); yValues.push_back(PxToCoord_y(65));
- //xValues.push_back(PxToCoord_x(668)); yValues.push_back(PxToCoord_y(63));
+
xValues.push_back(PxToCoord_x(759)); yValues.push_back(PxToCoord_y(85));
- //xValues.push_back(PxToCoord_x(775)); yValues.push_back(PxToCoord_y(143));
+
xValues.push_back(PxToCoord_x(832)); yValues.push_back(PxToCoord_y(159));
- //xValues.push_back(PxToCoord_x(874)); yValues.push_back(PxToCoord_y(139));
- //xValues.push_back(PxToCoord_x(815)); yValues.push_back(PxToCoord_y(36));
- //xValues.push_back(PxToCoord_x(848)); yValues.push_back(PxToCoord_y(53));
+
xValues.push_back(PxToCoord_x(938)); yValues.push_back(PxToCoord_y(101));
- xValues.push_back(PxToCoord_x(1100)); yValues.push_back(PxToCoord_y(62));
+ xValues.push_back(PxToCoord_x(1100)); yValues.push_back(PxToCoord_y(62));*/
+ xValues.push_back(PxToCoord_x(1240)); yValues.push_back(PxToCoord_y(70)); //Nuclei für Christians Vortrag
xValues.push_back(PxToCoord_x(1229)); yValues.push_back(PxToCoord_y(20));
xValues.push_back(PxToCoord_x(1258)); yValues.push_back(PxToCoord_y(29));
xValues.push_back(PxToCoord_x(1340)); yValues.push_back(PxToCoord_y(57));
- xValues.push_back(PxToCoord_x(1346)); yValues.push_back(PxToCoord_y(77));
- xValues.push_back(PxToCoord_x(1354)); yValues.push_back(PxToCoord_y(103));
- //xValues.push_back(PxToCoord_x(1432)); yValues.push_back(PxToCoord_y(101));
+ //xValues.push_back(PxToCoord_x(1346)); yValues.push_back(PxToCoord_y(77));
+ /*xValues.push_back(PxToCoord_x(1354)); yValues.push_back(PxToCoord_y(103));
+
xValues.push_back(PxToCoord_x(1731)); yValues.push_back(PxToCoord_y(115));
- xValues.push_back(PxToCoord_x(1915)); yValues.push_back(PxToCoord_y(77));
-*/
+ xValues.push_back(PxToCoord_x(1915)); yValues.push_back(PxToCoord_y(77));*/
+ // Addional Nuclei
+ //xValues.push_back(PxToCoord_x(668)); yValues.push_back(PxToCoord_y(63));
+ //xValues.push_back(PxToCoord_x(775)); yValues.push_back(PxToCoord_y(143));
+ //xValues.push_back(PxToCoord_x(874)); yValues.push_back(PxToCoord_y(139));
+ //xValues.push_back(PxToCoord_x(815)); yValues.push_back(PxToCoord_y(36));
+ //xValues.push_back(PxToCoord_x(848)); yValues.push_back(PxToCoord_y(53));
+ //xValues.push_back(PxToCoord_x(1432)); yValues.push_back(PxToCoord_y(101));
/*
//Experiment 7
@@ -491,7 +497,7 @@ public:
xValues.push_back(PxToCoord_x3(594)); yValues.push_back(PxToCoord_y3(680));
*/
-
+/*
//Inflow S1
xValues.push_back(PxToCoord_x4(1090)); yValues.push_back(PxToCoord_y4(1858));
xValues.push_back(PxToCoord_x4(1102)); yValues.push_back(PxToCoord_y4(1885));
@@ -515,7 +521,7 @@ xValues.push_back(PxToCoord_x4(1134)); yValues.push_back(PxToCoord_y4(2010));
xValues.push_back(PxToCoord_x4(1184)); yValues.push_back(PxToCoord_y4(2021));
xValues.push_back(PxToCoord_x4(1199)); yValues.push_back(PxToCoord_y4(1998));
xValues.push_back(PxToCoord_x4(1283)); yValues.push_back(PxToCoord_y4(2029));
-
+*/
}
@@ -530,6 +536,19 @@ xValues.push_back(PxToCoord_x4(1283)); yValues.push_back(PxToCoord_y4(2029));
temp = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-radius)));
yc = yc*temp;
}
+ if(phase >= 0)
+ {
+ RF yf = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-0.3)));
+ if (phase == 0)
+ {
+
+ yc = yc*yf;
+ }
+ else
+ {
+ yc = yc*(1-yf);
+ }
+ }
y[0] = yc;
}
};
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
index 84f8961d8363034e1e6f3427c9d1914643b155e4..834972d7f0fc3e23da52593d11bf9260236dc09d 100644
--- a/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
+++ b/dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh
@@ -82,21 +82,21 @@ public:
RF R1 = -Rspeed/param.eps*
(param.reaction.eval(ch.u)+param.phys.curvatureAlpha*(ch.xi-xi3))*ch.basis.phi[i];
RF J1 = param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma1 * (ch.gradxi*ch.basis.gradphi[i][0]);
-
+
// integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
r.accumulate(ch.pfspace,i,factor*(
(ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+J1
- +(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i]
- //-pf1*(v*gradphi[i][0]) weak transport
+ //+(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i]
+ -ch.pf*(ns.v*ch.basis.gradphi[i][0]) //weak transport
-R1
));
r.accumulate(ch.pf2space,i,factor*(
(ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma2 * (ch.gradxi2*ch.basis.gradphi[i][0])
- +(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i]
- //-pf2*(v*gradphi[i][0])
+ //+(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i]
+ -ch.pf2*(ns.v*ch.basis.gradphi[i][0])
));
r.accumulate(ch.xispace,i,factor*(
diff --git a/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh b/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a8c06c5fa0d62d03b8778fca76c9a23cf492884b
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh
@@ -0,0 +1,149 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+ class FEM_FFS_CHNS_R_ThinPlateFlow :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_FFS_CHNS_R_ThinPlateFlow<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FFS_CHNS_R_ThinPlateFlow<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+ {
+ private:
+ typedef typename CHContainer::field_type RF;
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CHLocalBasisCache<CHGFS> chCache;
+ VLocalBasisCache<NSGFS> vCache;
+ PLocalBasisCache<NSGFS> pCache;
+
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chLocal;
+ //mutable CHLView chOldLocal;
+
+ typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+ mutable NSLView nsOldLocal;
+
+ Param& param; // parameter functions
+
+ //FOR DETECTING ENTITY changes
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
+ mutable ENT OldEnt;
+
+ public:
+
+
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+
+
+
+ FEM_FFS_CHNS_R_ThinPlateFlow (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
+ NSGFS& nsGfs_, NSContainer& nsVOld_) :
+ chLocal(chGfs_,chV_),
+ //chOldLocal(chGfs_,chVOld_),
+ nsOldLocal(nsGfs_,nsVOld_),
+ param(param_)
+ {
+ }
+
+
+ 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
+ {
+
+ //Quadrature Rule
+ const int dim = EG::Entity::dimension;
+ const int order = 5; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chLocal.DoEvaluation(eg);
+ //chOldLocal.DoEvaluation(eg);
+ nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.lfs,Dune::Indices::_0));
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHVariables_Threephase_Extension<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+ NSVariables_P_VOld_Extension<RF,dim,LFSV> ns(vCache, pCache, ip.position(), lfsv, x.base(), nsOldLocal.lv, S);
+
+
+ RF pf3=1-ch.pf-ch.pf2;
+ RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
+ Dune::FieldVector<RF,dim> gradpf_f(0.0);
+ gradpf_f.axpy((1-2*param.delta),ch.gradpf+ch.gradpf2);
+
+ RF rho_f=param.phys.rho_f(ch.pf,ch.pf2);
+ RF rho_f_reg = rho_f + (param.phys.rho_1+param.phys.rho_2)*param.delta;
+
+ for (size_t i=0; i<ns.pspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // integrate: (div v)*phi
+ //r.accumulate(ns.pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor); Weak form
+ r.accumulate(ns.pspace,i, -1.0 * (pf_f*ns.div_v + (gradpf_f*ns.v) )* ns.pbasis.phi[i] * factor);
+ }
+ }
+ for (size_t i=0; i<ns.vfemspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // integrate: ( (rho+rhoOld)/2 * v_n+1 - rhoOld * v_n)/dt * phi + nu * Dv_n+1 : D phi - p_n+1 div phi
+
+ //std::cout << (param.nu * (jacv[k]*gradvphi[i][0]) - p * gradvphi[i][0][k]
+ //- param.rho*g[k]*vphi[i] )*factor << std::endl;
+
+ //RF lambda = 3;
+ //const auto v_nabla_v = v * jacv[k];
+ RF gradxi3k = -param.tw.Sigma3*(ch.gradxi[k]/param.tw.Sigma1 + ch.gradxi2[k]/param.tw.Sigma2);
+ RF S = 0;
+ S += - ch.xi2*(ch.gradpf[k]*pf_f - ch.pf * gradpf_f[k])/(pf_f*pf_f); //TODO!!!!!!! Other pf_f scaling!!!!!
+ S += - ch.xi *(ch.gradpf2[k]*pf_f - ch.pf2 * gradpf_f[k])/(pf_f*pf_f);
+ S += - 2*param.delta*pf3*(gradxi3k-ch.gradxi2[k]-ch.gradxi[k]);
+
+ //Stronger form - flow in interfaces possible
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k] - 2*param.delta*pf3*gradxi3k;
+
+ //Works but is not the one in the theory (same as above, w/o regularization)
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k];
+
+ S = S*param.SurfaceTensionScaling*ns.vbasis.phi[i];
+
+ r.accumulate(child(ns.vspace,k),i, ( rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt //TimeEvolution
+ //- ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k]) //Pressure in weak form
+ + pf_f * (ns.gradp[k] * ns.vbasis.phi[i]) //Pressure in strong form
+ + param.phys.eval_mu(pf_f) * (ns.jacv[k]*ns.vbasis.gradphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
+ + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i] //Velocity Dissipation in Solid
+ //- rho*g[k]*vphi[i] //Gravitation
+ - S //Surface Tension
+ + param.phys.eval_mu(pf_f) * 8 / (param.CellHeight * param.CellHeight) * ns.v[k]*ns.vbasis.phi[i]
+ )*factor);
+ if(k==0)
+ {
+ r.accumulate(child(ns.vspace,k),i, ( pf_f * (param.phys.outerMomentumXForce * ns.vbasis.phi[i]) // Momentum Forcing Term
+ )*factor);
+ }
+ } //for k
+ } // for i
+ } // for ip
+ } //alpha volume
+
+ };
diff --git a/src/paper_thinstrip_mixed/driver_mixed.hh b/src/paper_thinstrip_mixed/driver_mixed.hh
index 14342303e801d62e9f62f377040f3229b61826cc..7669eaf44c02aa27f0764000fdb38e8dc0ba39bb 100644
--- a/src/paper_thinstrip_mixed/driver_mixed.hh
+++ b/src/paper_thinstrip_mixed/driver_mixed.hh
@@ -12,6 +12,7 @@
#include "tmixed_chns_r_cahnhilliard.hh"
#include "tmixed_chns_navierstokes.hh"
#include "tmixed_upscaledp.hh"
+#include "tmixed_upscaledc.hh"
#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
@@ -202,8 +203,8 @@ void evaluateDxP(XGV& xgv, PGrad_DGF& pGradDgf, SimS& simS)
//_________________________________________________________________________________________________________________________
-template<typename GV, typename ES, typename TGV, typename T_FEM, typename V_FEM, typename P_FEM, typename PHI_FEM>
-void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
+template<typename GV, typename ES, typename TGV, typename T_FEM, typename V_FEM, typename P_FEM, typename PHI_FEM, typename C_FEM>
+void driver_mixed(GV& xgv, ES& x_es, const P_FEM& pFem, const C_FEM& cFem, const int xsize,
GV& ygv, ES& y_es, TGV& tgv, const T_FEM& tFem, const V_FEM& vFem, const PHI_FEM& phiFem,
Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
{
@@ -243,13 +244,11 @@ void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
//_________________________________________________________________________________________________________________________
// Make macro problems
- using XES = Dune::PDELab::NonOverlappingEntitySet<GV>;
- XES x_es(xgv);
-
typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ typedef Dune::PDELab::NoConstraints NCon;
ConstraintsAssembler conP();
- typedef Dune::PDELab::GridFunctionSpace<XES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
P_GFS pGfs(x_es,pFem);
auto blambda = [](const auto& x){return true;};
@@ -296,6 +295,47 @@ void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
pNewton.setLineSearchMaxIterations(10); // limit line search
pNewton.setAbsoluteLimit(1e-10);
+
+
+
+
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, C_FEM,Dune::PDELab::NoConstraints,VectorBackend> C_GFS;
+ C_GFS cGfs(x_es,cFem);
+ typedef Dune::PDELab::Backend::Vector<C_GFS, RF> C_V;
+ C_V cV(cGfs);
+ C_V cVOld(cGfs);
+ auto c0lambda = [&](const auto& xg)
+ {return 0.5+0.4*std::sin(4*3.14152*xg[0]);};
+ auto c0 = Dune::PDELab::
+ makeGridFunctionFromCallable(xgv,c0lambda);
+ Dune::PDELab::interpolate(c0,cGfs,cV);
+ cVOld = cV;
+ //typedef P_ParamAdaptor<RF, YSimS> P_PAdaptor;
+ //P_PAdaptor pPAdaptor(ySimS);
+ typedef TMIXED_UpscaledC<Parameters,P_PAdaptor,C_GFS,C_V,RF> C_LOP;
+ C_LOP cLop(param,pPAdaptor,cGfs,cVOld);
+ typedef Dune::PDELab::EmptyTransformation NoCC;
+ typedef Dune::PDELab::GridOperator<C_GFS, C_GFS, C_LOP, MBE,
+ RF, RF, RF, NoCC, NoCC> C_GO;
+ C_GO cGo(cGfs, cGfs, cLop, mbe);
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<C_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> C_LS;
+ C_LS cLs(cGo,500,100,3,1);
+ typedef Dune::PDELab::Newton<C_GO,C_LS,C_V> C_Newton;
+ C_Newton cNewton(cGo,cV,cLs);
+ cNewton.setReassembleThreshold(0.0); // always reassemble J
+ cNewton.setVerbosityLevel(3); // be verbose
+ cNewton.setReduction(1e-10); // total reduction
+ cNewton.setMinLinearReduction(1e-4); // min. red. in lin. solve
+ cNewton.setMaxIterations(25); // limit number of its
+ cNewton.setLineSearchMaxIterations(10); // limit line search
+ cNewton.setAbsoluteLimit(1e-10);
+
+
+
+
+
+
//_________________________________________________________________________________________________________________________
// prepare VTK writer
std::string filename=ptree.get("output.filename","output");
@@ -329,15 +369,23 @@ void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
P_DGF pDgf(pGfs,pV);
typedef Dune::PDELab::DiscreteGridFunctionGradient<P_GFS,P_V> PGrad_DGF;
PGrad_DGF pGradDgf(pGfs,pV);
+ typedef Dune::PDELab::DiscreteGridFunction<C_GFS,C_V> C_DGF;
+ C_DGF cDgf(cGfs,cV);
typedef Dune::PDELab::VTKGridFunctionAdapter<P_DGF> P_VTKF;
+ typedef Dune::PDELab::VTKGridFunctionAdapter<C_DGF> C_VTKF;
x_vtkwriter.addVertexData(std::shared_ptr<P_VTKF>(
new P_VTKF(pDgf,"p")));
+ x_vtkwriter.addVertexData(std::shared_ptr<C_VTKF>(
+ new C_VTKF(cDgf,"c")));
x_vtkwriter.write(param.time.t,Dune::VTK::base64);
//_________________________________________________________________________________________________________________________
+ bool FirstRun = true;
while (!timestepManager.isFinished())
{
+if(FirstRun)
+{
for(int i=1; i<xsize; i++)
{
ySimS[i].chVUpwind = ySimS[i-1].chV;
@@ -356,6 +404,7 @@ void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
}
ySimS[0].nsVUpwind = ySimS[xsize-1].nsV; //Periodic
+
//For Debug
//timestepManager.applyTimestep(vtkwriter);
//x_vtkwriter.write(param.time.t);
@@ -368,7 +417,11 @@ void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
{
ySimS[i].applyCh();
}
-
+FirstRun = false;
+}
+ std::cout << "Begin newtonapply Concentration:" << std::endl;
+ cNewton.apply();
+ cVOld =cV;
timestepManager.applyTimestep(vtkwriter);
x_vtkwriter.write(param.time.t);
}
diff --git a/src/paper_thinstrip_mixed/tmixed_equation.cc b/src/paper_thinstrip_mixed/tmixed_equation.cc
index 2ad9de3c5ef97318463f2212a557a9b4c721a5cb..ec4c9e5eae146a39ea4d9e61f207d42938ac8d7f 100644
--- a/src/paper_thinstrip_mixed/tmixed_equation.cc
+++ b/src/paper_thinstrip_mixed/tmixed_equation.cc
@@ -99,19 +99,20 @@ int main(int argc, char** argv)
using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
ES y_es(ygv);
+ ES x_es(xgv);
- const int k = 2;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> V_FEM;
- const int l = 2;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> PHI_FEM;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> P_FEM;
+ typedef Dune::PDELab::P0LocalFiniteElementMap<DF,double,1> C_FEM;
typedef Dune::PDELab::PkLocalFiniteElementMap<TGV,DF,double,1> Tensor_FEM;
V_FEM vFem(y_es);
- P_FEM pFem(y_es);
+ P_FEM pFem(x_es);
+ C_FEM cFem(Dune::GeometryTypes::cube(1));
PHI_FEM phiFem(y_es);
Tensor_FEM tensorFem(tgv);
- driver_mixed(xgv, pFem, A, ygv, y_es, tgv, tensorFem,vFem,phiFem, ptree, helper);
+ driver_mixed(xgv, x_es, pFem, cFem, A, ygv, y_es, tgv, tensorFem,vFem,phiFem, ptree, helper);
}
}
catch (Dune::Exception &e){
diff --git a/src/paper_thinstrip_mixed/tmixed_upscaledc.hh b/src/paper_thinstrip_mixed/tmixed_upscaledc.hh
new file mode 100644
index 0000000000000000000000000000000000000000..35ad4c3dad008f98d97cb1ed5e7abc9869230bf1
--- /dev/null
+++ b/src/paper_thinstrip_mixed/tmixed_upscaledc.hh
@@ -0,0 +1,120 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+
+
+
+
+template<typename Param, typename PAdaptor, typename CGFS, typename CContainer, typename RF>
+ class TMIXED_UpscaledC :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::FullSkeletonPattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<TMIXED_UpscaledC<Param, PAdaptor, CGFS, CContainer, RF> >,
+ public Dune::PDELab::NumericalJacobianSkeleton<TMIXED_UpscaledC<Param, PAdaptor, CGFS, CContainer, RF> >
+ {
+ private:
+ Param& param; // parameter functions
+ PAdaptor& pAdapt;
+
+ //typedef typename CGFS::Traits::Type FEM;
+ typedef typename CGFS::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
+ Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+
+ typedef Dune::PDELab::LocalFunctionSpace<CGFS> CLFS;
+ mutable CLFS cLfs;
+ typedef Dune::PDELab::LFSIndexCache<CLFS> CLFSCache;
+ mutable CLFSCache cLfsCache;
+ typedef typename CContainer::template ConstLocalView<CLFSCache> CView;
+ mutable CView cView;
+ mutable std::vector<typename LocalBasis::Traits::RangeFieldType> cOldLocal;
+
+ //FOR DETECTING ENTITY changes
+ typedef typename CGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
+ mutable ENT OldEnt;
+
+ public:
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+ enum { doPatternSkeleton = true};
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+ enum { doAlphaSkeleton = true};
+
+
+ TMIXED_UpscaledC (Param& param_, PAdaptor& pAdapt_, CGFS& cGfs_, CContainer& cVOld_) :
+ param(param_), pAdapt(pAdapt_),
+ cLfs(cGfs_), cLfsCache(cLfs), cView(cVOld_), cOldLocal(cGfs_.maxLocalSize())
+ {
+ }
+
+
+ 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
+ {
+ std::cout << "hi alpha" << std::endl;
+ auto c = x(lfsu,0);
+ if(eg.entity() != OldEnt)
+ {
+ OldEnt = eg.entity();
+ cLfs.bind(eg.entity());
+ cLfsCache.update();
+ cView.bind(cLfsCache);
+ cView.read(cOldLocal);
+ cView.unbind();
+ }
+ auto cOld = cOldLocal[lfsu.localIndex(0)];
+
+ r.accumulate(lfsv,0, (c-cOld)/param.time.dt*eg.geometry().volume());
+ std::cout << "bye alpha" << std::endl;
+ } //alpha volume
+
+
+ //! residual contribution from skeleton terms
+ template<typename IG, typename LFSU, typename X,
+ typename LFSV, typename R>
+ void alpha_skeleton (const IG& ig,
+ const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
+ const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
+ R& r_i, R& r_o) const
+ {
+ std::cout << "hi skeleton" << std::endl;
+ // inside and outside cells
+ auto cell_inside = ig.inside();
+ auto cell_outside = ig.outside();
+
+ // inside and outside cell geometries
+ auto insidegeo = cell_inside.geometry();
+ auto outsidegeo = cell_outside.geometry();
+
+ // cell centers in global coordinates
+ auto inside_global = insidegeo.center();
+ auto outside_global = outsidegeo.center();
+
+ // distance between the two cell centers
+ inside_global -= outside_global;
+ auto distance = inside_global.two_norm();
+
+ // face volume for integration
+ //auto facegeo = ig.geometry();
+ //auto face_volume = facegeo.volume();
+
+ // contribution to residual on inside and outside elements
+ auto dudn = (x_o(lfsu_i,0)-x_i(lfsu_o,0))/distance;
+ auto DiffusionCoeff = 0.5*(pAdapt.evalKf(inside_global) + pAdapt.evalKf(outside_global));
+ r_i.accumulate(lfsv_i,0,-dudn*DiffusionCoeff);
+ r_o.accumulate(lfsv_o,0, dudn*DiffusionCoeff);
+ std::cout << "bye skeleton" << std::endl;
+ }
+
+
+ };
diff --git a/src/summerschool/3p_driver.hh b/src/summerschool/3p_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..35dc6a9e611b9d5540e68040da43c71ee03264b3
--- /dev/null
+++ b/src/summerschool/3p_driver.hh
@@ -0,0 +1,261 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+#include<fstream>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+#include<dune/pdelab/function/callableadapter.hh>
+#include <dune/pdelab/common/functionutilities.hh>
+
+#include <dune/phasefield/localoperator/fem_ffs_chns_r_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_ffs_chns_r_thinplateflow.hh>
+#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/ffs_chns_r.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/boundary_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+#include "summerschool_timestepmanager.hh"
+#include "boundary_summerschool.hh"
+
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_ffs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+ typedef typename Grid::LeafGridView GV;
+ //std::bitset<2> periodic_directions(std::string("01"));
+ //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+ GV gv = grid.leafGridView();
+ const int dim = GV::dimension;
+ typedef double RF;
+
+
+ typedef Params_summerschool_3p< RF,
+ TripleWell<RF, DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> > >,//DoubleWell_poly<RF>,
+ Mobility_Degen<RF>,
+ Reaction_Calcite_Precipitation<RF>,
+ VelDissipation_Quadratic<RF> > Parameters;
+ Parameters param(ptree);
+ SummerschoolTimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+
+//_________________________________________________________________________________________________________________________
+// Make grid function spaces
+
+
+ typedef GFS_ffs_chns_r<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
+ GFS gfs(es,vFem,pFem,phiFem,phiFem);
+ typedef typename GFS::NS NS_GFS;
+ typedef typename GFS::CH CH_GFS;
+
+
+ //_________________________________________________________________________________________________________________________
+ //Constraints
+
+ typedef Bdry_ffs_chns_r_compositeV_2d<InflowOutflow_v_ss,
+ Dune::PDELab::DirichletConstraintsParameters,
+ InflowOutflow_p_ss,
+ InflowOutflow_phi_ss,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ InflowOutflow_phi_ss,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ InflowOutflow_u_ss> 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(gfs.ch);
+ NS_V nsV(gfs.ns);
+
+ //_________________________________________________________________________________________________________________________
+ //Initial Values
+
+ typedef Ini_ffs_chns_r< VelocityInitial_Inflow_ss<GV,RF,Parameters,dim>,
+ PInitial_Inflow<GV,RF,Parameters>,
+ PhiInitial_SummerSchoolPaper<GV,RF,Parameters,0>,
+ ZeroInitial<GV,RF,Parameters>,
+ PhiInitial_SummerSchoolPaper<GV,RF,Parameters,1>,
+ ZeroInitial<GV,RF,Parameters>,
+ UInitial_Inflow<GV,RF,Parameters> > Initial;
+ Initial initial(gv,param);
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+ //____________________________________________________________________________________________________________________
+ // Refine Initial Data
+
+ for(int i=0; i<ptree.get("Refinement.maxRefineLevel",2); i++)
+ {
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_3P_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+ RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+ rInd.calculateIndicator(gv,chV,false);
+ rInd.markGrid(grid,2);
+
+ // do refinement
+ auto chT = transferSolutions(gfs.ch,4,chV);
+ auto nsT = transferSolutions(gfs.ns,4,nsV);
+ Dune::PDELab::adaptGrid(grid, chT,nsT);
+
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+ }
+
+
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+
+ CH_V chVOld = chV;
+ NS_V nsVOld = nsV;
+
+ //_________________________________________________________________________________________________________________________
+ // LocalOperators
+
+ typedef FEM_FFS_CHNS_R_ThinPlateFlow<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+ NS_LOP nsLop(param, gfs.ch, chV, gfs.ns, nsVOld);
+
+ typedef FEM_FFS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
+ CH_LOP chLop(param, 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().
+
+
+ //_________________________________________________________________________________________________________________________
+ // Construct discrete Grid Functions
+ typedef DGF_ffs_chns_r<GFS, NS_V, CH_V> DGF;
+ DGF dgf(gfs, nsV, 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,1);
+ Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
+ dgf.addToVTKWriter(vtkwriter);
+ vtkwriter.write(param.time.t,Dune::VTK::base64);
+
+
+
+
+
+ //_________________________________________________________________________________________________________________________
+ printVector(chV,10,"chV");
+ printVector(nsV,10,"nsV");
+
+
+ //____________________________________________________________________________________________________________________
+ //int i_it=0;
+
+
+ while (!timestepManager.isFinished())
+ {
+
+ // Grid Operators
+ 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,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 ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+ <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,1000,200,3,1);
+
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+ <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,4000,400,3,1);
+
+ // Nonlinear solver
+ typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
+ PDESolver1 ns_newton(nsGo,nsV,ls_ns);
+ ns_newton.setParameters(ptree.sub("NS_Newton"));
+
+ typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
+ PDESolver2 ch_newton(chGo,chV,ls_ch);
+ ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
+
+ //_________________________________________________________________________________________________________________________
+ // do time step
+
+
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
+
+ int result;
+ result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+
+ if(result <0)
+ {
+ return;
+ }
+ else if(result == 0)
+ {
+ nsV = nsVOld;
+ chV = chVOld;
+ continue;
+ }
+
+ //_________________________________________________________________________________________________________________________
+
+
+ chVOld = chV;
+ nsVOld = nsV;
+
+ //_________________________________________________________________________________________________________________________
+
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_3P_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+ RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+ rInd.calculateIndicator(gv,chV,false);
+ rInd.markGrid(grid,2);
+
+ printVector(chV,3,"chV");
+ printVector(nsV,3,"nsV");
+
+ // do refinement
+ std::cout << "adapt grid and solution... ";
+ auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
+ auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
+ Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
+ std::cout << " ...adapt ended" << std::endl;
+
+ printVector(chV,3,"chV");
+ printVector(nsV,3,"nsV");
+
+ //vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+
+ // recompute constraints
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+
+
+ //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
+ }
+}
diff --git a/src/summerschool/CMakeLists.txt b/src/summerschool/CMakeLists.txt
index c08acc1adbb08518a80770c55d707e9785a50045..6f94a3a0209dcf0775ca7828c2030476e31c8391 100644
--- a/src/summerschool/CMakeLists.txt
+++ b/src/summerschool/CMakeLists.txt
@@ -1,2 +1,3 @@
add_executable(summerschool summerschool.cc)
+#add_executable(summerschool3p summerschool3p.cc)
dune_symlink_to_source_files(FILES grids ss_calcite.ini)
diff --git a/src/summerschool/boundary_summerschool.hh b/src/summerschool/boundary_summerschool.hh
index 88fadea977e044c28693e9c069332c4303aa5282..4cb1e146e917a3da8ca40a86d0fafa7c96d51dc6 100644
--- a/src/summerschool/boundary_summerschool.hh
+++ b/src/summerschool/boundary_summerschool.hh
@@ -6,7 +6,7 @@
#include <dune/phasefield/initial_fn.hh>
#define GRID_EPS_SS 1e-7
-#define GRID_WIDTH_SS 7 // 1.5 for summer school //2 for contactpointtest //6 for full summer school, 7 for v2
+#define GRID_WIDTH_SS 1.5 // 1.5 for summer school //2 for contactpointtest //6 for full summer school, 7 for v2
#define GRID_BOTTOM_SS 0 //0 For summer school //0.35 for contactpointtest //0.35 for hainesjump
#define GRID_TOP_SS 0.125 //0.125 for summer school //1 for contactpointtest //0.65 for hainesjump
diff --git a/src/summerschool/ss_calcite.ini b/src/summerschool/ss_calcite.ini
index 98622da5b7cffcaa91aa1d1256233241ed769a4d..ab02eda0da5bdcce1bd667e18784561f39e46d89 100644
--- a/src/summerschool/ss_calcite.ini
+++ b/src/summerschool/ss_calcite.ini
@@ -1,7 +1,8 @@
[Physics]
-uAst = 27.1 #Molar density of CaCO3, in mol/l
+uAst = 0.271 #Molar density of CaCO3, in mol/l #27.1, for Christian Presentation: 0.271
uEquilibrium = 0.00014 #Solubility of CaCO3 in water, in mol/l
-ReactionRate = 10
+ReactionRate = 10 #For Reaction
+#ReactionRate = 0 #For No Reaction (Threephase)
BulkRate = 0.00005716#0.00001141 #0.00005716 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s) # Was 0.00001796
D = 0.000804 #of CO3^2- in Water see Zeebe,2011, in mm^2/s
rho1 = 0.997 #Water, in 10^-6 kg/mm^3
@@ -9,31 +10,38 @@ rho2 = 2.71 #Kalcite, in 10^-6 kg/mm^3
mu=1.01 #Viscosity of Water, in mPa*s = 10^-6 kg/mm
VelDissipationRate = 1000000
CurvatureAlpha=0
+#SurfaceTensionScaling=100 #72000 #SurfaceTension in 10^-6 kg/s^2 #Only for Threephase
[domain]
-filename = grids/summerschool_full_v2.msh
+#filename = grids/summerschool_full_v2.msh
+filename = grids/summerschool.msh
CellHeight=0.085
[output]
-filename = finalresult
-
+filename = precip_christian_v2
+#filename = tpexperimental
[Phasefield]
-eps=0.002
-delta=0.01
-Mobility=0.00002#0.001 * 0.1 #New *0.2
+eps=0.002 #0.002
+delta=0.03
+Mobility=0.00002#0.001 * 0.1 #New *0.2 For 2 phase
+#Mobility = 0.02 # For 3phase
+#MobilityMin = 0.002 # For 3phase
DoubleWellDelta=0.08
DoubleWellGamma=0.06
+#Sigma1=1 #Only for Threephase
+#Sigma2=1 #Only for Threephase
+#Sigma3=10 #Only for Threephase
[Time]
tMax = 100000 #25000 * 4
-dt = 0.0015
-dtMax = 0.1
-dtMin = 0.0001
-saveintervallstart = 100
-saveintervall = 119.99
-reactionOnlyDt = 20 #NEW
-reactionOnlyStart = 20 #NEW
+dt = 0.0015 #0.0015
+dtMax = 0.05 #0.1
+dtMin = 0.0001 #0.0001
+saveintervallstart = 40
+saveintervall = 2
+reactionOnlyDt = 2 #NEW
+reactionOnlyStart = 4 #NEW
#dtIncreaseFactor=1.2
#dtDecreaseFactor=0.5
@@ -54,9 +62,10 @@ x3Center=0.8
y3Center=0.225
[Refinement]
-etaRefine = 0.30
+etaRefine = 0.33
etaCoarsen = 0.15
-maxRefineLevel = 7
+#maxRefineLevel = 7
+maxRefineLevel = 6
minRefineLevel = 0
[NS_Newton]
diff --git a/src/summerschool/ss_calcite_backup.ini b/src/summerschool/ss_calcite_backup.ini
new file mode 100644
index 0000000000000000000000000000000000000000..0b4dbee6fc066e6d0d1444efce4d6e99041d9145
--- /dev/null
+++ b/src/summerschool/ss_calcite_backup.ini
@@ -0,0 +1,83 @@
+[Physics]
+uAst = 27.1 #Molar density of CaCO3, in mol/l
+uEquilibrium = 0.00014 #Solubility of CaCO3 in water, in mol/l
+#ReactionRate = 10 #For Reaction
+ReactionRate = 0 #For No Reaction (Threephase)
+BulkRate = 0.00005716#0.00001141 #0.00005716 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s) # Was 0.00001796
+D = 0.000804 #of CO3^2- in Water see Zeebe,2011, in mm^2/s
+rho1 = 0.997 #Water, in 10^-6 kg/mm^3
+rho2 = 2.71 #Kalcite, in 10^-6 kg/mm^3
+mu=1.01 #Viscosity of Water, in mPa*s = 10^-6 kg/mm
+VelDissipationRate = 1000000
+CurvatureAlpha=0
+SurfaceTensionScaling=100 #72000 #SurfaceTension in 10^-6 kg/s^2 #Only for Threephase
+
+[domain]
+#filename = grids/summerschool_full_v2.msh
+filename = grids/summerschool.msh
+CellHeight=0.085
+
+[output]
+#filename = threephase
+filename = tpexperimental
+
+[Phasefield]
+eps=0.002 #0.002
+delta=0.03
+#Mobility=0.00002#0.001 * 0.1 #New *0.2 For 2 phase
+Mobility = 0.02 # For 3phase
+MobilityMin = 0.002 # For 3phase
+DoubleWellDelta=0.08
+DoubleWellGamma=0.06
+Sigma1=1 #Only for Threephase
+Sigma2=1 #Only for Threephase
+Sigma3=10 #Only for Threephase
+
+[Time]
+tMax = 100000 #25000 * 4
+dt = 0.0005 #0.0015
+dtMax = 0.001 #0.1
+dtMin = 0.00005 #0.0001
+saveintervallstart = 100
+saveintervall = 119.99
+reactionOnlyDt = 20 #NEW
+reactionOnlyStart = 20 #NEW
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+MeanInflow=1.882 #mm/s (in mean-> *1.5 for max velocity)
+uConst=0.000315003 #FROM INFLOW SIMULATION
+uInflow=0.000315003 #FROM INFLOW SIMULATION
+Radius=0.02
+
+Radius1=0.015
+x1Center=0.8
+y1Center=-0.1
+Radius2=0.015
+x2Center=0.6
+y2Center=0.11
+Radius3=0.015
+x3Center=0.8
+y3Center=0.225
+
+[Refinement]
+etaRefine = 0.33
+etaCoarsen = 0.15
+#maxRefineLevel = 7
+maxRefineLevel = 6
+minRefineLevel = 0
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 1
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 1
diff --git a/src/summerschool/summerschool3p.cc b/src/summerschool/summerschool3p.cc
new file mode 100644
index 0000000000000000000000000000000000000000..162b5c50064b9780072a937d2e4cd0f2e4d103e7
--- /dev/null
+++ b/src/summerschool/summerschool3p.cc
@@ -0,0 +1,99 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include<iostream>
+// dune includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.hh>
+//#if HAVE_UG
+//#include<dune/grid/uggrid.hh>
+//#endif
+#if HAVE_DUNE_ALUGRID
+#include<dune/alugrid/grid.hh>
+#include<dune/alugrid/dgf.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+// pdelab includes
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+
+#include "3p_driver.hh"
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+
+
+
+int main(int argc, char** argv)
+{
+ #if !HAVE_SUPERLU
+ std::cout << "Error: These examples work only if SuperLU is available." << std::endl;
+ exit(1);
+ #endif
+
+ #if HAVE_UG
+ try{
+ // Maybe initialize Mpi
+ Dune::MPIHelper&
+ helper = Dune::MPIHelper::instance(argc, argv);
+ if(Dune::MPIHelper::isFake)
+ std::cout<< "This is a sequential program." << std::endl;
+ else
+ std::cout << "Parallel code run on "
+ << helper.size() << " process(es)" << std::endl;
+
+ // open ini file
+ Dune::ParameterTree ptree;
+ Dune::ParameterTreeParser ptreeparser;
+ ptreeparser.readINITree("ss_calcite.ini",ptree);
+ ptreeparser.readOptions(argc,argv,ptree);
+
+ typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+ //typedef Dune::UGGrid<2> Grid;
+ typedef Grid::ctype DF;
+
+ std::vector<int> boundary_index_map;
+ std::vector<int> element_index_map;
+
+ std::string filename = ptree.get("domain.filename",
+ "turbtube2d.msh");
+ Dune::GridFactory<Grid> factory;
+ if (helper.rank()==0)
+ Dune::GmshReader<Grid>::read(factory,filename,boundary_index_map,element_index_map,true,false);
+ std::shared_ptr<Grid> gridp(factory.createGrid());
+
+ typedef Grid::LeafGridView GV;
+ GV gv = gridp->leafGridView();
+
+ std::cout << "OK" << std::endl;
+
+ using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+ ES es(gv);
+
+ const int k = 2;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k-1> P_FEM;
+ const int l = 2;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+ V_FEM vFem(es);
+ P_FEM pFem(es);
+ PHI_FEM phiFem(es);
+
+ driver_ffs(*gridp, es,vFem,pFem,phiFem, ptree, helper);
+
+
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 1;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ return 1;
+ }
+ #endif
+}