diff --git a/src/paper_thinstrip_full/3p_driver_monolithic.hh b/src/paper_thinstrip_full/3p_driver_monolithic.hh index bfb11a84f90787beaf3ae9e25c234091f274970c..3c2e1c28636c71943f4cd227b38f4eb9483d7753 100644 --- a/src/paper_thinstrip_full/3p_driver_monolithic.hh +++ b/src/paper_thinstrip_full/3p_driver_monolithic.hh @@ -37,7 +37,7 @@ void driver_3p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem, TripleWell<RF,DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> > >, Mobility_Const<RF>, Reaction_Thinstrip<RF>, - VelDissipation_Quadratic<RF> > Parameters; + VelDissipation_Quadratic_Shifted<RF> > Parameters; Parameters param(ptree); TimestepManager timestepManager(ptree.sub("Time"),param.time,2); diff --git a/src/paper_thinstrip_full/chns.ini b/src/paper_thinstrip_full/chns.ini index 83b64852f10c947f1edcf6bb1d2c95eb907b303e..bcab18aa9cde7def748616cf18426941f7e0f3a7 100644 --- a/src/paper_thinstrip_full/chns.ini +++ b/src/paper_thinstrip_full/chns.ini @@ -8,6 +8,7 @@ ReactionRate=0.5 #ReactionLimiter = 1000 SurfaceTensionScaling=0.1 VelDissipationRate = 1000 +#VelDissipationMax = 0.9 CurvatureAlpha=0.001 [domain] @@ -31,7 +32,7 @@ Sigma2=1 Sigma3=1 [Time] -tMax = 10 +tMax = 1 dt = 0.0015 dtMax = 0.01 dtMin = 0.00015 diff --git a/src/paper_thinstrip_full/thinstrip_boundary.hh b/src/paper_thinstrip_full/thinstrip_boundary.hh index 86f9af6307617909669ce599038403aa01acdfe2..34b36a27f91bdffbf1401dc929345c3346961e47 100644 --- a/src/paper_thinstrip_full/thinstrip_boundary.hh +++ b/src/paper_thinstrip_full/thinstrip_boundary.hh @@ -102,7 +102,7 @@ public: RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[1] - (dSolid + d1+d_sine)*pref.YScaling)); if(Phase == 1) { - y = rf2; + y = rf*rf2; return; } y=rf*(1-rf2); diff --git a/src/paper_thinstrip_mixed/chns.ini b/src/paper_thinstrip_mixed/chns.ini index 82fa6b1240b53db9e10d771d85b87a99cd0aa755..79e3111d6b973a6750a21dc76ac3cfd907af7369 100644 --- a/src/paper_thinstrip_mixed/chns.ini +++ b/src/paper_thinstrip_mixed/chns.ini @@ -4,11 +4,12 @@ D = 0.02 rho1 = 1.0 rho2 = 1.0 mu=0.005 -ReactionRate=1.5 +ReactionRate=0.5 #ReactionLimiter = 1000 SurfaceTensionScaling=0.1 -VelDissipationRate = 100 -CurvatureAlpha=0#0.001 +VelDissipationRate = 1000 +#VelDissipationMax = 0.9 +CurvatureAlpha=0.001 #0 [domain] level = 0 @@ -19,21 +20,21 @@ YResolution = 200 XResolution = 100 [output] -filename = output2 +filename = comparison_hyperbolic [Phasefield] eps=0.03 -delta=0.00001#0.03 +delta= 0.03 #0.00001 DoubleWellDelta=0.05 -Mobility=1#0.01 +Mobility=0.04 #1 #0.01 Sigma1=1 Sigma2=1 Sigma3=1 [Time] -tMax = 100 -dt = 0.0001 -dtMax = 0.0001 +tMax = 1 +dt = 0.001 +dtMax = 0.001 dtMin = 0.000015 saveintervall = 0.005 #dtIncreaseFactor=1.2 @@ -41,7 +42,7 @@ saveintervall = 0.005 [Initial] # FlowVyy= 10# pConst/(Domain_width * mu) -pConst=0.1 +pConst=0.05 uConst=0.5 uInflow=0.3 diff --git a/src/paper_thinstrip_mixed/driver_mixed.hh b/src/paper_thinstrip_mixed/driver_mixed.hh index e1ad831826dd2d91eaebc00d639495fcab61e28c..57fb59ae8d79f6fabc00f9a17f0600c6ab8674cd 100644 --- a/src/paper_thinstrip_mixed/driver_mixed.hh +++ b/src/paper_thinstrip_mixed/driver_mixed.hh @@ -174,6 +174,7 @@ public: Dune::FieldVector<double,1> integral; Dune::PDELab::integrateGridFunction(productFunction,integral,2); kf = integral[0]; + std::cout << "Kf = " << kf << std::endl; } void updatePhiC() @@ -283,7 +284,7 @@ void evaluateDxP(XGV& xgv, PGrad_DGF& pGradDgf, SimS& simS) for(int i=0; i<xs; i++) { pGradDgf.evaluate(*it,it->geometry().local(it->geometry().center()), pgrad ); - pgrad[0] = pgrad[0] - 1; //TODO make more general with pressure drop/distance + pgrad[0] = pgrad[0] - 0.05; //TODO make more general with pressure drop/distance simS[i].dxp = pgrad[0]; simS[i].dxpUpwind = pgradUpwind[0]; it++; @@ -324,8 +325,8 @@ void driver_mixed(XGV& xgv, X_ES& x_es, const P_FEM& pFem, const C_FEM& cFem, co TripleWell<RF,DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> > >, Mobility_Const<RF>, Reaction_Thinstrip<RF>, - VelDissipation_Quadratic<RF> - //VelDissipation_Quadratic_Shifted<RF> + //VelDissipation_Quadratic<RF> + VelDissipation_Quadratic_Shifted<RF> > Parameters; Parameters param(ptree); TimestepManager timestepManager(ptree.sub("Time"),param.time,2); diff --git a/src/paper_thinstrip_mixed/tmixed_boundary.hh b/src/paper_thinstrip_mixed/tmixed_boundary.hh index 7da76192b58a59063023d31327606d436a9570c1..60af9cd2f3aa4e0841e51de1f9c1fd8e8304f45a 100644 --- a/src/paper_thinstrip_mixed/tmixed_boundary.hh +++ b/src/paper_thinstrip_mixed/tmixed_boundary.hh @@ -168,17 +168,17 @@ public: { eps = params.eps; dSolid = params.initial.ptree.get("dSolid",0.3); - d1 = params.initial.ptree.get("d1",0.3); + d1 = params.initial.ptree.get("d1",0.4); } inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const { //bool inMiddle = (pref.TempXpos > 0.4) && (pref.TempXpos < 0.6); //RF mod = 0.1*std::sin(2*3.141592*pref.TempXpos);//inMiddle?0.3:0; - RF mod = 0; + RF d_sine = 0.15*sin(pref.TempXpos*2*M_PI); //std::cout << "TempXPos" << pref.TempXpos; - RF rf = pref.tw.dw.eval_shape(1/eps*(x[0] - dSolid -mod)); - RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[0] - dSolid - d1 )); + RF rf = pref.tw.dw.eval_shape(1/eps*(x[0] - dSolid)); + RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[0] - (dSolid + d1 + d_sine))); if(Phase == 1) { y = rf*rf2; diff --git a/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh b/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh index 7c2c3a701cda39d6d4d9a2ddad23fe1d5f9deb6a..93a2e2dc29776b6a680ca9a8545966bdb69c2e49 100644 --- a/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh +++ b/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh @@ -119,9 +119,9 @@ public: //RF Rspeed = std::max(param.tw.eval_q(ch.pf,ch.pf2,pf3) - 0.1, 0.0); RF Rspeed = param.tw.eval_q(ch.pf,ch.pf2,pf3); - auto cMod = c -0.2*std::sin(2*3.14152*xpos); + //auto cMod = c -0.2*std::sin(2*3.14152*xpos); RF R1 = -Rspeed/param.eps* - (param.reaction.eval(cMod)+param.phys.curvatureAlpha*(ch.xi-xi3))*ch.basis.phi[i]; + (param.reaction.eval(c)+param.phys.curvatureAlpha*(ch.xi-xi3))*ch.basis.phi[i]; // RF R1 = 0; diff --git a/src/paper_thinstrip_mixed/tmixed_upscaledp.hh b/src/paper_thinstrip_mixed/tmixed_upscaledp.hh index a40e910ec1bb000865f1c54c70afd2773b8f783c..12bb241326ab689adb4e7588bee454c9bc049f89 100644 --- a/src/paper_thinstrip_mixed/tmixed_upscaledp.hh +++ b/src/paper_thinstrip_mixed/tmixed_upscaledp.hh @@ -75,7 +75,7 @@ template<typename Param, typename FEM, typename RF> p += x(lfsv,i)*phi[i]; } auto xg = eg.geometry().global(ip.position()); - gradp[0] = gradp[0] - 1; //TODO: Make more general for other pressure jumps/distances + gradp[0] = gradp[0] - 0.05; //TODO: Make more general for other pressure jumps/distances //auto xg = eg.geometry().center(); //std::cout << "pos: " << xg[0] << " kf: " << kfresult-kfold << std::endl; auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());