diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index e928570a090ff858cf701ea748d422f6030ed955..0fdafb5b37a2ffb75f21cfd423353bd5f69a64c7 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -64,6 +64,7 @@ public:
lastsave = time.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)
@@ -87,6 +88,30 @@ public:
return 1;
}
+
+ // RETURN VALUE: -1 -> ABORT, 0-> REDO TIMESTEP, 1-> OK
+ template<typename NS_NEWTON, typename CH_NEWTON, typename VTKWRITER>
+ int doTimestepNoAccept(NS_NEWTON& ns_newton, CH_NEWTON& ch_newton, VTKWRITER& vtkwriter)
+ {
+ try
+ {
+ if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
+ ns_newton.apply();
+
+ if (verb > 1) std::cout << "= Begin newtonapply CH:" << std::endl;
+ ch_newton.apply();
+
+ }
+ catch(...)
+ {
+ return rejectTimestep();
+ }
+
+ vtkwriter.write(time.t-1e-5,Dune::VTK::base64);
+
+ return 1;
+ }
+
// RETURN VALUE: -1 -> ABORT, 0-> REDO TIMESTEP, 1-> OK
template<typename NS_NEWTON, typename VTKWRITER>
int doOnlyNSTimestep(NS_NEWTON& ns_newton, VTKWRITER& vtkwriter)
@@ -165,7 +190,7 @@ public:
if( time.t < saveintervallstart || time.t>lastsave+saveintervall )
{
lastsave = time.t;
- vtkwriter.write(time.t,Dune::VTK::ascii);
+ vtkwriter.write(time.t,Dune::VTK::base64);
}
}
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index 88bb72ff5f64d26a4ab729369a68c0c698e0ea14..601bf0548a0ddea8fc068544e3df1d7a043fe32f 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -400,8 +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)-3.5;
+ return 0.5+(px_x - 150)/((RF)299);
+ //return 0.5+(px_x - 150)/((RF)299)-3.5; //-> Christians Vortrag
}
inline static RF PxToCoord_y(const RF px_y)
{
@@ -446,7 +446,9 @@ public:
//Experiment 1
- xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75));
+ xValues.push_back(PxToCoord_x(114)); yValues.push_back(PxToCoord_y(81)); //Extra Nucleus, to check the influence on Nucleus 1
+
+ xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75)); //Nucleus 1
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));
@@ -461,7 +463,7 @@ public:
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(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));
diff --git a/src/paper_thinstrip_full/3p_driver.hh b/src/paper_thinstrip_full/3p_driver.hh
index ac38b762193160262f0752e0e8b5ca48f87ec839..eabac9e8049c427b414b9396984d454ea88265a2 100644
--- a/src/paper_thinstrip_full/3p_driver.hh
+++ b/src/paper_thinstrip_full/3p_driver.hh
@@ -56,16 +56,16 @@ void driver_3p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem,
//_________________________________________________________________________________________________________________________
//Constraints
typedef DomainSize<RF> DS;
- DS dS(0,param.YScaling,0,2,1e-8);
+ DS dS(0,1,0,ptree.get("domain.len",1),1e-8);
typedef Bdry_ffs_chns_r_compositeV_2d_Explicit_Domain<
BottomBoundary<DS>, // v_x
AllBoundary<DS>, // v_y
- LeftRightBoundary<DS>, // p
- LeftBoundary<DS>, // phi1
+ NoBoundary<DS>, // p
+ NoBoundary<DS>, // phi1
NoBoundary<DS>, // mu1
- LeftBoundary<DS>, // phi2
+ NoBoundary<DS>, // phi2
NoBoundary<DS>, // mu2
- LeftBoundary<DS>> // u
+ NoBoundary<DS>> // u
Bdry;
Bdry bdry(dS);
gfs.updateConstraintsContainer(bdry);
@@ -189,7 +189,19 @@ void driver_3p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem,
interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
- int result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+ int result = timestepManager.doTimestepNoAccept(ns_newton,ch_newton,vtkwriter);
+ if(result <0)
+ {
+ return;
+ }
+ else if(result == 0)
+ {
+ nsV = nsVOld;
+ chV = chVOld;
+ continue;
+ }
+
+ result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
//int result = timestepManager.doOnlyNSTimestep(ns_newton,vtkwriter);
if(result <0)
{
diff --git a/src/paper_thinstrip_full/chns.ini b/src/paper_thinstrip_full/chns.ini
index 75f37777cc21f4bb2b55d81e85bc289a8601a60b..c0f2d28d0006d263c465cc477f99d23e68fb4285 100644
--- a/src/paper_thinstrip_full/chns.ini
+++ b/src/paper_thinstrip_full/chns.ini
@@ -12,12 +12,14 @@ CurvatureAlpha=0.001
[domain]
level = 0
-filename = grids/rectangle_periodic_coarse.msh
-YScalingFactor = 0.5
+filename = grids/thinstrip_periodic_2.msh
+len = 2
+periodic = 01
+#YScalingFactor = 0.5
YScaling = 0.5
[output]
-filename = output2
+filename = output4
[Phasefield]
eps=0.03
@@ -29,7 +31,7 @@ Sigma2=1
Sigma3=1
[Time]
-tMax = 30
+tMax = 300
dt = 0.015
dtMax = 0.15
dtMin = 0.0015
diff --git a/src/paper_thinstrip_full/thinstrip_boundary.hh b/src/paper_thinstrip_full/thinstrip_boundary.hh
index e7f453d3d149bc5a8d64ccf7253ce0a156bd93a9..e241f05b0c31516fa8f5f2cc443f0361bdbf4f34 100644
--- a/src/paper_thinstrip_full/thinstrip_boundary.hh
+++ b/src/paper_thinstrip_full/thinstrip_boundary.hh
@@ -133,8 +133,8 @@ public:
Dune::PDELab::AnalyticGridFunctionBase<Traits,VelocityInitial_InOutFlow<GV,RF,Params,dim> > (gv)
, time(params.time)
{
+ //TODO: Use right scaling (now small coordinates)
v = params.initial.ptree.get("pConst",(RF)0.01)/(params.YScaling*params.YScaling) /2 /*Domain Width*/ /params.phys.mu;
- yS = params.YScaling;
}
template <typename T>
@@ -144,7 +144,7 @@ public:
inline void evaluateGlobal(const DomainType & x, RangeType & y) const
{
- y[0] = std::max(0.0, v/2.0*(x[1]-0.3*yS)*((2*GRID_TOP_THINSTRIP-0.3)*yS-x[1]));
+ y[0] = std::max(0.0, v/2.0*(x[1]-0.3)*((2*GRID_TOP_THINSTRIP-0.3)-x[1]));
}
};
@@ -154,7 +154,6 @@ class PInitial_InOutFlow :
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PInitial_InOutFlow<GV,RF,Params> >
{
RF pConst;
- RF yS;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
@@ -162,16 +161,15 @@ public:
Dune::PDELab::AnalyticGridFunctionBase<Traits,PInitial_InOutFlow<GV,RF,Params> > (gv)
{
pConst = params.initial.ptree.get("pConst",(RF)0.01);
- yS = params.YScaling;
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
y = 0;
- if(x[0]<1e-7)
+ /*if(x[0]<1e-7)
{
y=pConst/(yS*yS);
- }
+ }*/
}
};
@@ -181,7 +179,6 @@ class PhiInitial_Thinstrip_Layers :
Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Thinstrip_Layers<GV,RF,Params,Phase> >
{
RF eps;
- RF yS;
RF dSolid;
RF d1;
const Params & pref;
@@ -196,13 +193,12 @@ public:
eps = params.eps;
dSolid = params.initial.ptree.get("dSolid",0.3);
d1 = params.initial.ptree.get("d1",0.3);
- yS = params.YScaling;
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
- RF rf = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid*yS));
- RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid*yS - d1*yS));
+ RF rf = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid));
+ RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid - d1));
if(Phase == 1)
{
y = rf2;
diff --git a/src/paper_thinstrip_full/thinstrip_equation.cc b/src/paper_thinstrip_full/thinstrip_equation.cc
index f89e7c46644e94811823601cf60b78b25ecd741c..a80993f5eec3039cf062c1f08190c86919264ac3 100644
--- a/src/paper_thinstrip_full/thinstrip_equation.cc
+++ b/src/paper_thinstrip_full/thinstrip_equation.cc
@@ -25,7 +25,7 @@
// Main program with grid setup
//===============================================================
-template<typename RF>
+/*template<typename RF>
class ExampleThinstripTransform
: public Dune :: AnalyticalCoordFunction< RF, 2, 2, ExampleThinstripTransform<RF> >
{
@@ -84,7 +84,7 @@ class StrechingThinstripTransform
y[0] = (x[0]-c)*(S-2*q)/(1-2*q)+c;
}
}
-};
+};*/
int main(int argc, char** argv)
@@ -129,14 +129,26 @@ int main(int argc, char** argv)
std::string filename = ptree.get("domain.filename",
"turbtube2d.msh");
+// Unstructured Grid
+ Dune::GridFactory<Grid> factory;
+ Dune::GmshReader<Grid>::read(factory,filename,true,true);
+ Dune::GridFactory<Grid>::WorldMatrix matrix;
+ std::cout << "a" << std::endl;
+ matrix[0][0] = 1; matrix[0][1] = 0;
+ matrix[1][0] = 0; matrix[1][1] = 1;
+ std::cout << "b" << std::endl;
+ Dune::GridFactory<Grid>::WorldVector shift;
+ shift[0] = ptree.get("domain.len",1);
+ shift[1] = 0;
+ std::cout << "c" << std::endl;
+ factory.insertFaceTransformation(matrix,shift );
+ std::shared_ptr<Grid> gridp(factory.createGrid());
+ std::cout << "d" << std::endl;
+
+
- //Dune::GridFactory<Grid> factory;
- //Dune::GmshReader<Grid>::read(factory,filename,true,true);
- //std::shared_ptr<Grid> gridp(factory.createGrid());
- double A = ptree.get("domain.OneOverYGridWidth", 2);
- double B = ptree.get("domain.YScalingFactor",1.0);
// Transforming big rectangle by shortening y-Direction
// Dune::FieldVector<DF,2> lowerLeft(0);
@@ -150,44 +162,50 @@ int main(int argc, char** argv)
//ExampleThinstripTransform<double>* tfp = new ExampleThinstripTransform<double>(ptree.get("domain.YScalingFactor",1.0));
// Transforming small rectangle by elongating in x-Direction
- Dune::FieldVector<DF,2> lowerLeft(0);
- Dune::FieldVector<DF,2> upperRight(0);
- upperRight[0]=B;
- upperRight[1]=B/A;
- std::array<unsigned int,2> elements;
- elements[0]=ptree.get("domain.CellsInYDirection", 32)*A;
- elements[1]=ptree.get("domain.CellsInYDirection", 32);
+ // double A = ptree.get("domain.OneOverYGridWidth", 2);
+ // double B = ptree.get("domain.YScalingFactor",1.0);
+ // Dune::FieldVector<DF,2> lowerLeft(0);
+ // Dune::FieldVector<DF,2> upperRight(0);
+ // upperRight[0]=B;
+ // upperRight[1]=B/A;
+ // std::array<unsigned int,2> elements;
+ // elements[0]=ptree.get("domain.CellsInYDirection", 32)*A;
+ // elements[1]=ptree.get("domain.CellsInYDirection", 32);
+ //
+ // typedef Dune::GeometryGrid<Grid,StrechingThinstripTransform<double>> GGrid;
+ // StrechingThinstripTransform<double>* tfp = new StrechingThinstripTransform<double>(1/B);
+// End of possible Transformations
- typedef Dune::GeometryGrid<Grid,StrechingThinstripTransform<double>> GGrid;
- StrechingThinstripTransform<double>* tfp = new StrechingThinstripTransform<double>(1/B);
- std::shared_ptr<Grid> gridp(Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, elements));
- std::shared_ptr<GGrid> ggridp(new GGrid(gridp.get(),tfp));
+// Structured Grid for Transformations
+ // std::shared_ptr<Grid> gridp(Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, elements));
+ // std::shared_ptr<GGrid> ggridp(new GGrid(gridp.get(),tfp));
+ // typedef GGrid::LeafGridView LGV;
- typedef GGrid::LeafGridView LGV;
+ typedef Grid::LeafGridView LGV;
std::bitset<2> periodic_directions(ptree.get("domain.periodic","00"));
typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
- GV gv(ggridp->leafGridView(),periodic_directions);
+ GV gv(gridp->leafGridView(),periodic_directions);
//typedef Grid::LeafGridView GV;
//GV gv = gridp->leafGridView();
std::cout << "OK" << std::endl;
- using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+ using ES = Dune::PDELab::AllEntitySet<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;
+ const int l = 1;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
V_FEM vFem(es);
P_FEM pFem(es);
PHI_FEM phiFem(es);
- driver_3p(*ggridp, gv, es,vFem,pFem,phiFem, ptree, helper);
+ driver_3p(*gridp, gv, es,vFem,pFem,phiFem, ptree, helper);
}
}
catch (Dune::Exception &e){
diff --git a/src/paper_thinstrip_full/thinstrip_full_chns_navierstokes.hh b/src/paper_thinstrip_full/thinstrip_full_chns_navierstokes.hh
index ba0bc4b116bd14f8d3f1ccfbde5367a8e196b489..3773f6a8e04b442348e2f491bcf4db9f544751c8 100644
--- a/src/paper_thinstrip_full/thinstrip_full_chns_navierstokes.hh
+++ b/src/paper_thinstrip_full/thinstrip_full_chns_navierstokes.hh
@@ -88,6 +88,8 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
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);
+ Dune::FieldVector<RF,dim> outerGradP(0.0);
+ outerGradP[0]=-param.gradP;
gradpf_f.axpy((1-2*param.delta),ch.gradpf+ch.gradpf2);
RF rho_f=param.phys.rho_f(ch.pf,ch.pf2);
@@ -129,7 +131,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
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
+ + pf_f * ((ns.gradp[k]+outerGradP[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
diff --git a/src/paper_thinstrip_full/thinstrip_parameters.hh b/src/paper_thinstrip_full/thinstrip_parameters.hh
index 1a820a024cb506b6e49dcf2e37229b9a5c564d3c..b444058c94f18de16d1ad58b19e32d94deef4854 100644
--- a/src/paper_thinstrip_full/thinstrip_parameters.hh
+++ b/src/paper_thinstrip_full/thinstrip_parameters.hh
@@ -9,16 +9,20 @@ class Params_thinstrip : public Params_ffs_r<Number, TW, Mobility, Reaction, Vel
{
public:
const Number YScaling;
+ const Number len;
+ const Number gradP;
Params_thinstrip(const Dune::ParameterTree& ptree) :
Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>(ptree),
- YScaling(ptree.get("domain.YScaling",1.0))
+ YScaling(ptree.get("domain.YScaling",1.0)),
+ len(ptree.get("domain.len",1.0)),
+ gradP(ptree.get("Initial.pConst",1.0) / len) //TODO: Right scaling!
{
ModifyParameters();
std::cout << "YScaling " << YScaling << " Eps " << this->eps << std::endl;
}
- void ModifyParameters()
+ void ModifyParameters() //Now TODO, we are in small coordinates!
{
this->eps *= YScaling;
this->mobility.mobility *= YScaling*YScaling;
diff --git a/src/rectangle_periodic.geo b/src/rectangle_periodic.geo
index 434fe3bd1fa1b93bd64daa17977a405cfc98a55c..d2957b01e3b9db5d99e1d71f44f1011d776ada2b 100644
--- a/src/rectangle_periodic.geo
+++ b/src/rectangle_periodic.geo
@@ -1,10 +1,11 @@
// mesh width associated with points
-lc = 0.015;
+lc = 0.03;
+len = 2;
Point(1) = {0, 0, 0, lc};
-Point(2) = {0.9375, 0, 0, lc};
-Point(3) = {0.9375, 0.75, 0, lc};
-Point(4) = {0, 0.75, 0, lc};
+Point(2) = {len, 0, 0, lc};
+Point(3) = {len, 1, 0, lc};
+Point(4) = {0, 1, 0, lc};
Line(1) = {2,1};
Line(2) = {3,2};
diff --git a/src/summerschool/ss_calcite.ini b/src/summerschool/ss_calcite.ini
index 3a57355e0a40f28b2604ec2a0034367282b168dd..3665b2bebedec33a44d7dc1246a21727ec07a04d 100644
--- a/src/summerschool/ss_calcite.ini
+++ b/src/summerschool/ss_calcite.ini
@@ -13,14 +13,14 @@ 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
+filename = grids/summerschool_full_v2.msh
+#filename = grids/summerschool.msh
CellHeight=0.085
[output]
#filename = precip_christian_v2
#filename = tpexperimental
-filename = flowonly
+filename = extra_nucleus
[Phasefield]
eps=0.002 #0.002