//
// Phase Field - Implicit + characteristics
//
// The problem:
//
// $u_t + (u\cdot\nabla)u - \nabla\cdot(\nu(\phi)Du) + \nabla p = f + \theta g$ 
// $\nabla \cdot u = 0$
// $\theta_t + u\cdot\nabla\theta - \nabla\cdot(\kappa(\phi)\nabla\theta) 
//                                                            = \nu(\phi)Du:Du$
// $\phi_t + u\cdot\nabla\phi - a\Delta\phi = -bH(\phi) + c|\nabla\phi|\theta$
// 
// in $\Omega\times(0,T)$
//
// $u = u_\Gamma$
// $\kappa(\phi)\partial_n\theta = h$
// $\partial_n\phi = 0$
//
// on $\partial\Omega\times(0,T)$

// $u|_{t=0} = u^0$
// $\theta|_{t=0} = \theta^0$
// $\phi|_{t=0} = \phi^0$
//
// in $\Omega$
//
// The functions:
//
// $\nu(s) = \nu_0 + \nu_1 (1 - s)$
// $H(s) := s(1-s)(1-2s)$)
// $\kappa(s) = \kappa_0 + \kappa_1 (1 - s)$
//
// The variables:
//
// $u, p$: velocity field and pressure
// $\phi$: phase field function
// $\theta$: temperature
//
//
// At each time step, we solve the following stationary problem in 
// $[t^n, t^n + \Delta t]$:
//
// - ${1 \over \Delta t} (u^{n + 1} - u^n \otimes X^n)
//       -\nabla\cdot(\nu(\phi^n)Du^{n + 1}) + \nabla p^{n + 1} 
//       = f^{n + 1} + \theta^n g$
//   $\nabla \cdot u^{n + 1} = 0$
//   together with boundary conditions.
//   
// - ${1 \over \Delta t} (\theta^{n + 1} - \theta^n \otimes X^n)
//       -\nabla\cdot(\kappa(\phi^n)\nabla\theta^{n + 1}) 
//       = \nu(\phi^n)Du^{n + 1}:Du^{n + 1}$
//   together with boundary conditions.
//   
// - ${1 \over \Delta t} (\phi^{n + 1} - \phi^n \otimes X^n)
//       -a\Delta\phi^{n + 1})
//       = -bH(\phi^{n}) + |\nabla\phi^{n}|\theta^{n + 1}$
//   together with boundary conditions.
//
//
//
//
//

// Mesh 1
border aa(t=0., 2*pi) {    x = 5.*cos(t);    y = 5.*sin(t);  };
border cc(t=0., 2*pi) {    x = -3.+0.2*cos(t);    y = 0.2*sin(t);  };
mesh Th = buildmesh(aa(30)+cc(12));
plot(Th,wait=1);
int ireg = Th(-3.,0.).region;

// Mesh 2
// real aa = 101, bb = 102;
// border a1(t=0,10) {x=t; y=0;label=aa;}
// border a2(t=0,10) {x=10; y=t;label=aa;}
// border a3(t=0,10) {x=10-t; y=10;label=aa;}
// border a4(t=0,10) {x=0; y=10-t;label=aa;}
// border b1(t=0,2) {x=1+t; y=4; label=bb;}
// border b2(t=0,2) {x=3; y=4+t; label=bb;}
// border b3(t=0,2) {x=3-t; y=6; label=bb;}
// border b4(t=0,2) {x=1; y= 6-t; label=bb;}
// int n=2;
// mesh Th = buildmesh(a1(5)+a2(5)+a3(5)+a4(5)+b1(2)+b2(2)+b3(2)+b4(2));
// plot(Th,wait=1);
// iregcont = Th(2.,5.).region;

// Constants and functions to furnish
real nu00 = 1.e-1, nu01 = 1.e1;        // data for the viscosity
real kappa00 = 1.e-1, kappa01 = 1.e1;  // data for the heat diffusion
real a=1.e-1, b=1., c=1.;              // constants in the equation for $\phi$
real alpha00 = 15., uinf00 = 1.;       // data for the boundary velocity
real u00 = 0., v00 = 0.;               // data for the initial velocity
real theta00=0.;                       // data for the initial temperature
real phi0a=1., phi0b=0.;               // data for the initial phase field
real fext00 = 0., gext00 = 0.;         // data for the external force
real g100 = 1.e-1, g200 = 1.e-1;       // data for the caloric force
real h00 = 5.e-1;                      // data for the boundary heat flux
real T=10.;                             // final time
int ntsteps = 500;                     // number of time steps

// Auxiliar constants and variables
real dnu0nu1 = 2.*(nu00 + nu01), dnu01 = 2.*nu01,
		al00 = alpha00*pi/180., uGamma = uinf00*cos(al00), vGamma = uinf00*sin(al00),
		ka0ka1 = kappa00 + kappa01, dt = T/ntsteps, alphaB = 1./dt;
int n;
real epsm = 1., epsphi = 1.e-5;

func u0f = u00;                                  // initial $u_1$
func v0f = v00;                                  // initial $u_2$
func fextf = fext00;                             // $f_1$
func gextf = gext00;                             // $f_2$
func g1f = g100;                                 // $g_1$
func g2f = g200;                                 // $g_1$
//func hf = h00*(x >= 0.);                                   // $h$
func hf = 0.;                                    // $h$
func theta0 = theta00;
// func fextf = -fext0*(Th(x,y).region == iregcont)*y;    // $f_1$
// func gextf = gext0*(Th(x,y).region == iregcont)*(x+3); // $f_2$
func phi0 = phi0b + (phi0a-phi0b) * (Th(x,y).region == ireg);

// Spaces, physical constants and parameters
fespace Vh(Th,P1b); // P1-bubble for the velocity
fespace Ph(Th,P1);  // P1 for the pressure
fespace Zh(Th,P2);  // P2 for the temperature and the phase field

// Functions
Vh u=u0f, v=v0f, uu, vv, uold, vold, f, g, fsecond, gsecond;
Ph p=0., pp, pold, psi=0., w, dnucoef, kacoef;
Zh theta=theta0, thetat, thetaold, zeta, thsecond;
Zh phi=phi0, phit, phiold, xi, eta, phsecond;

// The equations
//
// Motion
// Stationary problems

problem qstokes([u,v,p],[uu,vv,pp]) = int2d(Th)(alphaB*(u*uu+v*vv) 
        + dnucoef*(2.*dx(u)*dx(uu)+(dy(u)+dx(v))*(dy(uu)+dx(vv))+2.*dy(v)*dy(vv))
		+ dx(p)*uu+dy(p)*vv+pp*(dx(u)+dy(v))-1.e-10*p*pp)
		- int2d(Th)(fsecond*uu + gsecond*vv)
// Instead: - int2d(Th)((alphaB*uold-eta*(uold*dx(uold)+vold*dy(uold)))*uu
//                     +(alphaB*vold-eta*(uold*dx(vold)+vold*dy(vold)))*vv)
        + on(aa, u=uGamma, v=vGamma);

problem heat(theta,thetat) =
		int2d(Th)(alphaB*theta*thetat
		+ kacoef*(dx(theta)*dx(thetat)+dy(theta)*dy(thetat)))
        - int2d(Th)(thsecond*thetat)
        - int1d(Th,aa)(hf*thetat);


problem phase(phi,phit) = int2d(Th)(alphaB*phi*phit
		+ a*(dx(phi)*dx(phit)+dy(phi)*dy(phit)))
		- int2d(Th)(phsecond*phit);

macro tau(U)(N.x*uGamma-N.y*vGamma) //
problem pb4psi(psi,w,init=n,solver=LU) = int2d(Th)(dx(psi)*dx(w)+dy(psi)*dy(w))
		- int2d(Th)((dx(v) - dy(u))*w) 
		// - int1d(Th,aa)(tau(U)*w)
+ on(aa,psi = 0.);


// Initial plots
plot([u,v],fill=1,value=1,cmm="INITIAL VELOCITY FIELD",wait=true);
plot(theta,fill=1,value=1,cmm="INITIAL TEMPERATURE",wait=true);
plot(phi,fill=1,value=1,cmm="INITIAL PHASE FIELD",wait=true);

// Iterates
for(int n=0; n<ntsteps; n++){
	// Actualizing old functions
	uold = u; vold = v; pold=p;
	thetaold = theta;
	phiold = phi;
	// For material time derivatives
	f = convect([u,v],-dt,uold); g = convect([u,v],-dt,vold);
	zeta = convect([u,v],-dt,thetaold);
	xi = convect([u,v],-dt,phiold);
	
	// Solving the stationary problem for $u$ from t to t + dt	
	dnucoef = dnu0nu1 - dnu01*phi;
	fsecond = alphaB*f + fextf + theta*g1f;
	gsecond = alphaB*g + gextf + theta*g2f;
	qstokes;
	// plot(p,fill=1,value=1);
	// plot(u);
	// plot([u,v]);
	
	// Solving the stationary problem for $\theta$ from t to t + dt
	zeta = convect([u,v],-dt,thetaold);
	kacoef = ka0ka1 - kappa01*phi;
	thsecond = alphaB*zeta + 
			kacoef*(2.*dx(u)*dx(u)+(dy(u)+dx(v))*(dy(u)+dx(v))+2.*dy(v)*dy(v));
	heat;
	// plot(theta,fill=1,value=1);

	// Solving the stationary problem for $\phi$ from t to t + dt
	xi = convect([u,v],-dt,phiold);
	eta = dx(phi)*dx(phi)+dy(phi)*dy(phi);
	phsecond = alphaB*xi - b*phi*(1.-phi)*(1.-2.*phi) + c*sqrt(eta)*theta;
	phase;
	plot(phi,fill=1,value=1);
	
//	xi = convect([u,v],-dt,phiold);
//	xi = alphaB*xi;
//
//	while(epsm >= epsdir){
//		if(m > mmax) break;
//		phim = phi; 
//		phsecond = xi-b*phi*(1.-phi)*(1.-2.*phi)
//				+ c*sqrt{dx(phi)*dx(phi)+dy(phi)*dy(phi)}*theta;
//		phase;
//		epsm = error(phi - phim);
//	};
//	plot(phi,fill=1);
	
	// Mesh adaptation:
	if((n+1) % 20 == 0){
		Th = adaptmesh(Th,phi,nbvx=500);
		plot(Th, wait=true);
		plot([u,v],wait=true);
	}
};

// Final figures:
// plot(u,wait=1,value=1,ps="PhFu-I1.eps");
plot([u,v],wait=1,ps="PhFvel-I1.eps");
plot(p,wait=1,fill=1,value=1,ps="PhFp-I1.eps");
plot(theta,wait=1,fill=1,value=1,ps="PhFtheta-I1.eps");
plot(phi,wait=1,fill=1,value=1,ps="PhFphi-I1.eps");