// file CH-G-OS.edp
//
// OPTIMAL CONTROL FOR THE HEAT PDE - GRADIENT WITH OPTIMAL STEP (AND PROJECTION)
//

verbosity=0;

//
// NUMERICAL VALUES TO PROVIDE
//
real T=10.;
real kappa=1., y00=1.;
int C0=100, C1=99, C2=98; // the labels
int Ntimes=20;

real umin=0., umax=100.;
real acoef=.95, one=1., yd0=1.; // good values of acoef: <= .95
real rhofix=1.;
real errGO=1.e-5, errmesh=1.e-1;
int itermax=500;

//
// AUXILIAR VARIABLES
//
int kiter=0, jiter=0, nmesh=1, itermm1=itermax-1, Nt1=Ntimes+1;

real tau=T/Ntimes, xnu=kappa*tau, bcoef=1.-acoef, abcoef=acoef/bcoef,
    abc=rhofix*one, error, errn, aux1, aux2, rho;
real[int] Uerror(1000);

//
// DESCRIPTION OF THE BOUNDARIES
//
// outer boundary ($\Gamma_0$):
border C00(t=0,2*pi){x=5*cos(t); y=5*sin(t); label=C0;}

// inner boundary ($\Gamma_1$):
border C11(t=0.,1.){ x=2.-t;  y=3.;      label=C1;}
border C12(t=0.,1.){ x=1.;    y=3.-6.*t;  label=C1;}
border C13(t=0.,1.){ x=1.+t;  y=-3.;     label=C1;}
border C14(t=0.,1.){ x=2.;    y=-3.+6.*t; label=C1;}

//
// FIGURES (I)
//
// boundaries:
plot (                C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20),
                    wait=true);

//
// MESH
//
mesh Th=buildmesh(    C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20));

//
// FIGURES (II)
//
// initial mesh:
plot(Th, cmm=" FIRST MESH", wait=true);

//
// THE FINITE ELEMENT SPACE AND FUNCTIONS
//
fespace Vh(Th,P1);
// fespace Vh(Th,P2);
// Vh u=0.,v,fff,f0=f00,u0=u00,uold,udif;
Vh ystate, pstate, zstate, sol, phi, yd;
Vh udif, fff;
Vh[int] un(Nt1), u(Nt1), des(Nt1);

for (int n=0;n<Ntimes;n++)
{
    un[n]=0.;
}

int iref0=Th(2.1,3.1).region;
int iref1=Th(1.9,2.9).region;

//
// SOME FUNCTIONS
//
func unom=one*(region==iref1);
func abom=abc*(region==iref1);
func ydfun=yd0*(x*x-y*y);
func y0=y00;

yd=ydfun;

// -------------------------------------
// THE SYSTEMS
// -------------------------------------
problem states(sol,phi,init=jiter) = int2d(Th)(sol*phi
        + xnu*(dx(sol)*dx(phi)+dy(sol)*dy(phi)))
        - int2d(Th)(fff*phi)
        + on(C0,sol=0.);
// -------------------------------------

// -------------------------------------
// THE FIRST PROJECTOR
// -------------------------------------
// gvec = max(fvec,umin)
//
func real[int] PUad(real[int] & fvec)
{
	real[int] gvec(fvec.n);
	for (int ifv=0;ifv<fvec.n;ifv++)
	{
		gvec[ifv] = max(fvec[ifv],umin);
	}
	
	return gvec;
}
// -------------------------------------

/*
// -------------------------------------
// THE SECOND PROJECTOR
// -------------------------------------
// gvec = max(min(fvec,umax),umin)
//
func real[int] PUad(real[int] & fvec)
{
	real[int] gvec(fvec.n);
	for (int ifv=0;ifv<fvec.n;ifv++)
	{
		gvec[ifv] = min(max(fvec[ifv],umin),umax);
	}
	
	return gvec;
}
// -------------------------------------
*/

//
// SOLVING THE PROBLEM
//
for (int n=0;n<Ntimes;n++)
{
    u[n]=un[n];
}


// ITERATES

for (int k=0;k<itermax;k++)
{
	kiter = k;
	
	ystate=y0;
	for (int n=1;n<=Ntimes;n++)
	{
		fff=ystate+tau*u[n];
		states;
		jiter=1;
		ystate=sol;
	    plot(ystate, cmm="STATE - ITERATE "+kiter, value=true,
	        fill=true, wait=false);
	}
	plot(ystate, cmm="STATE - ITERATE "+kiter, value=true,
	    fill=true, wait=true);
	    
	aux1=0.;	 
	fff = ystate-yd;
	for (int n=Ntimes;n>0;n--)
	{
		states;
		pstate=sol;
	    plot(pstate, cmm="ADJOINT STATE - ITERATE "+kiter, value=true,
	        fill=true, wait=false);
	    un[n]=u[n];
	    des[n]=abom*(acoef*pstate+bcoef*un[n]);
	    aux1=aux1+int2d(Th)(des[n]*des[n]);
//        u[n]=u[n]-abom*des[n];
//        u[n][]=PUad(u[n][]);
// 	      udif = u[n]-un[n];
//	      errn = udif[].l2;
//	      error=error+errn*errn;
		fff=pstate;
	}
	plot(pstate, cmm="ADJOINT STATE - ITERATE "+kiter, value=true,
	    fill=true, wait=true);

	zstate=0.;
	for (int n=1;n<=Ntimes;n++)
	{
		fff=zstate+tau*des[n];
		states;
		zstate=sol;
// 	      plot(ystate, cmm="STATE - ITERATE "+kiter, value=true,
//          fill=true, wait=false);
	}
    
    aux2=int2d(Th)(zstate*zstate);
    rho=aux1/(acoef*aux2+bcoef*aux1);

    error=0.;
	for (int n=1;n<=Ntimes;n++)
	{
		u[n]=u[n]-rho*des[n];
//      u[n][]=PUad(u[n][]);
 	    udif = u[n]-un[n];
	    errn = udif[].l2;
	    error=error+errn*errn;
        plot(u[n], cmm="CONTROL - ITERATE "+kiter, value=true,
	        fill=true, wait=false);
	}
    plot(u[Ntimes], cmm="CONTROL - ITERATE "+kiter, value=true,
	    fill=true, wait=true);
	    
	error=sqrt(error);
	Uerror[kiter]=error;
	cout << " ITER = " << kiter << " ERROR = " << error << endl;
	if (error <= errGO) break;
	
/*
	if ((kiter+1) % 100 == 0)
	{
		nmesh = nmesh+1;
		jiter = -1;
		Th = adaptmesh(Th,u,err=errmesh);
		plot(Th, cmm=" MESH NO. "+nmesh, wait=true);
		u=u;
	};
*/

	if (kiter == itermm1) cout << " ATTENTION: TOO MANY ITERATES" << endl;
	
	for (int n=1;n<=Ntimes;n++)
	{
        un[n]=u[n];
	}
	
	jiter = jiter+1;
}

//
// DISPLAY, RECORD AND PRINT
//

for (int n=0;n<Ntimes;n++)
{
    plot(u[n], cmm="CONTROL - ITERATE "+kiter, value=true,
	    fill=true, wait=false);
}
plot(u[Ntimes], cmm="CONTROL - ITERATE "+kiter, value=true,
	fill=true, wait=true);

plot(ystate, value=true, fill=true, cmm="STATE - FINAL ITERATE, FINAL TIME",
    wait=1);
plot(pstate, value=true, fill=true, 
    cmm="ADJOINT STATE - FINAL ITERATE, INITIAL TIME", wait=1);


// WRITING

/*
{ // save solution
	ofstream file("sol.txt");
	file << u[] << endl;
} // close the file (end block)
// { // read
// ifstream file("sol.txt");
// file >> u[] ;
// } // close reading file (end block)

kiter1 = kiter+1;
for (int k=0;k<kiter1;k++)
{
	cout << " FIXED-POINT ITERATE AND ERROR: " << k << " - "<< Uerror[k] << endl;
}

if (kiter > 1)
{
	aux = log(Uerror[kiter-1]/Uerror[kiter])/log(2.);
    cout << " FIXED-POINT RATE =             " << aux << endl;
}

//
// FIGURES (III)
//
// the solution:
plot(u, wait=true, value=true, fill=true);
// plot(u, wait=true, value=true, fill=true, ps="nonsymsol.eps");
*/