// model for mutation-advection-selection,  Chisholm, Lorenzi et al. 2015
// with NIPG
// IMEX scheme: semi-backward differentiation formula, SBDF
// domain
int n = 40; // size
mesh Th = square(n,n,flags=1);
Th = adaptmesh(Th, hmax=1./60.);
// FE space - discontinuous Galerkin
fespace Vh(Th,P2dc);
Vh u,v, uold,uold1;

// define operators: Gradient
macro Grad(u) [ dx(u),dy(u) ] //
// normal times adv. field
macro nn() (N.y*vv) //
macro uup(u) (mean(u)-(nn<0)*jump(u)/2) //

// model parameters
// drug dose
real m = 0.15,vvv=1.2e-2;
// function beta(x,y) - delta
func bb = .03 + .25*y+.05*(1-x);
func dd = m*(.15 + 1.3*(1-x)) + .0;

func vv = -vvv*m*(x<=.5); // adv. field

// parameters
real alpha = 1.e-5;
real I0,I1;
// integration
real dt = Th.hmin/(m*vvv), t=0.;
dt = fmin(.5,.1*dt);

// weak formulation 
// computing stiffness matrix & vector
//===================================================
// non-symmetric interior penalty Galerkin (NIPG)
// parameters for DG penalisation terms
real sigma0 = 1.7, epsilon = 1, b0 = 1;
// variational formulations
varf a1(u,v) = int2d(Th)(u*v/dt); // time derivative

// Laplace operator (Riviere, ch. 4)
varf a2(u,v)= int2d(Th)( alpha*(dx(u)*dx(v) + dy(u)*dy(v)) )
	- intalledges(Th) ( (1-nTonEdge)*sigma0/(lenEdge^b0) *jump(u)*jump(v) )
	+intalledges(Th)( (1-nTonEdge)*( alpha* mean( N.x*dx(u)+N.y*dy(u) ) )*jump(v) )
	- intalledges(Th) ( epsilon*((1-nTonEdge)*( alpha* mean( N.x*dx(v)
	+N.y*dy(v) )  )*jump(u)) ) ;

// advection operator (Riviere, ch. 4)
varf b(u,v) = int2d(Th)( (vv*dy(v) )*u ) 
	+ intalledges(Th) ( (1-nTonEdge) * nn *uup(u)*jump(v) ) 
	- int1d(Th) ( (nn>0)*nn*u*v );

// reaction terms
varf ll1(u,v) = int2d(Th) (bb *u*v) ;
varf ll2(u,v) = int2d(Th) (dd *u*v) ;
// build stiffness matrices: bilinear part
matrix A1 =a1(Vh,Vh),A2=a2(Vh,Vh);
// build load vectors: linear part
matrix L1 = ll1(Vh,Vh), L2 = ll2(Vh,Vh), BB = b(Vh,Vh);;

Vh rhs,rhs0 ,rhs1, err;

int iter =800;
real tol = 1.e-4, dtol=1;
// initial data
real x0=.2,y0=.8,c1=200.;
uold = c1*exp(-(x-x0)^2/(2*.01^2))*exp(-(y-y0)^2/(2*.01^2))*((x-x0)^2+(y-y0)^2<.1^2); uold=abs(uold); // u{1}
//
// need an auto-starting method for time-integration
// find u{2} by IMEX-theta-scheme
// compute integral
I0 = int2d(Th)(uold);
matrix A0 = A1+2./3.*A2;
set (A0,init=1,solver=UMFPACK);
rhs[] = A1*uold[];
rhs1[] = A2*uold[]; rhs1 = -1./3.*rhs1;
rhs[] += rhs1[];
rhs1[] = L1*uold[];
rhs1 = (1-I0)*rhs1;
rhs[] += rhs1[];
rhs1[]= L2*uold[]; 
rhs[] -=rhs1[];
rhs1[]= BB*uold[]; 
rhs[] +=rhs1[];
u[]=A0^-1*rhs[];

matrix A = 3./2.*A1+A2;
set (A,init=1,solver=UMFPACK);
// loop
for (int citer = 0; citer <= iter; citer++)
{
	t+=dt;
// initialise RHS of the variational problem
// compute integrals
	I0 = int2d(Th)(u);
	I1 = int2d(Th)(uold);
//	cout << "Int(u) = "<<I0<< endl;
// add summands
	rhs[]=A1*u[]; // 2*<u{n}, phi> + 2*dt*int u{n} <u{n}, phi> 
	rhs = 2.*rhs;
	rhs1[]= L1*u[]; rhs1 = 2.*(1.-I0)*rhs1; // 2*dt*<(b-d)u{n}, phi >
	rhs[] += rhs1[];
	rhs1[]=L2*u[]; rhs1 = 2.*rhs1;
	rhs[] -=rhs1[];
	rhs1[]= BB*u[]; rhs1 =2*rhs1;
	rhs[] +=rhs1[];
	rhs1[] = A1*uold[]; rhs1 = -0.5*rhs1;// -1/2*<u{n-1}, phi>+dt* int u{n-1}< u{n-1}, phi >
	rhs[] += rhs1[];
	rhs1[]=L1*uold[]; rhs1 = (1-I1)*rhs1; // -dt* <(b-d)*u{n-1}, phi>
	rhs[] -= rhs1[];
	rhs1[]= BB*uold[];
	rhs[] -=rhs1[];
	rhs1[]=L2*uold[]; // -dt* <(b-d)*u{n-1}, phi>
	rhs[] += rhs1[];
// solve explicit part
	uold = u;
	u[]=A^-1*rhs[];
	plot(u,value=1,fill=1,cmm=" t="+t );
//	err[] = u[] - uold[];
//	dtol = err[].linfty;
//	cout << "dtol = "<<dtol<< endl;
}