//
//    Stratigraphic model with diffusive sediment transport 
//    Modelize the formation of clastic sedimentary basins at large space and time scales 
//
//    Single lithology Model 
//
//
//  d_t h(x,t) + div(grad(psi(b(x,t)))) = 0 on (0,Lx)x(0,tf)
//
//  h(x,0) = hinit(x)   on (0,Lx)
//
//  grad(psi).n = g0 for x=0,  grad(psi).n = g1 for x = Lx  
//
//  b(x,t) = hmer(t) - h(x,t)
//
//
//  Finite Volume discretization using TPFA fluxes on unstructured meshes 
//
//  Ouputs: * h(x,t), b(x,t)
//
//   
//          * hs(x,t) = min_(t<=q<=tf) h(x,q) (sediment layers at each time t, taking into account erosions)
//


clear

//data
Lx=2;
//space discretization
N=100;
Nint = N-1;
Nbound = 2;
dx=Lx/N;

//time discretization
tf = 1.5; // final simulation time
dtinit = 0.02; // initial time step
dtmax = 0.02; // maximum time step 
Ndtmax = 10000; // maximum number of time step

//Newton convergence 
Newtmax = 10; // maximum number of Newton iterations 
eps=1.0E-6; //  stopping criteria 


Km = 1;
Kc = 10;

g0 = -20;
g1 = 0;

function s = f_hmer(t)
    s = 25 + 5*cos(12*t) ;
endfunction

function s = f_psi(u)
    if (u<0) then
        s = Kc*u;
    else 
        s = Km*u;
    end
endfunction

function s = f_psip(u)
    if (u<0) then
        s = Kc;
    else 
        s = Km;
    end
endfunction


function s = f_hinit(x)
    s = 25*exp(-8*x/Lx) + 10;
endfunction


function R = residual(h,h0,hmer)

for m=1:N
  //  R(m) = 
end

for i=1:Nint
 
end

for i=1:Nbound
   
end

endfunction


function A = Jacobian(h,hmer)

A=sparse([],[],[N,N]);

for m=1:N
  //  A(m,m) = 
end

for i=1:Nint

end


endfunction



//data structure for the uniform 1D mesh of size dx of the domain (0,Lx)
// cells m = 1:N
for m=1:N
volume(m) = dx
X(m) = (m-1/2)*dx
end  

//Interior faces: i = 1:N-1
for i=1:Nint

surface(i) = 1
cellint(i,1) = i 
cellint(i,2) = i+ 1 
end

// Boundary
cellbound(1) = 1;
cellbound(2) = N;
fbound(1) = g0;
fbound(2) = g1;


//transmissibilities of interior faces
for i=1:Nint
Tint(i) = surface(i)/abs(X(cellint(i,1))-X(cellint(i,2)))
end

/////////////////////////////////


//initial sediment thickness
for k=1:N
 h0(k)= f_hinit(X(k))
end

dt = dtinit;
time(1) = 0;
timestep(1) = dt;

scf(2);
title('h')
plot (X,h0,'r')  


//Initialization of Newton ------------------------
h = h0;

//time loop -----------------------------------------------------------
t = 0;
n = 1;
newt = 1;
timer()
while ((t<tf)&(n<Ndtmax))
n = n + 1;
t = t + dt
hmer = f_hmer(t);

////////////////////
// to be completed 
////////////////////

//Newton loop -----------------------------------
itn = 1;
xn(itn) = newt;
yn(itn) = 1;

while ((normR/normR0>eps)&(itn<Newtmax)) do 


////////////////////
// to be completed 
////////////////////


itn = itn+1
newt = newt + 1
xn(itn) = newt
yn(itn) = normR/normR0

end //end Newton loop ---------------------


if (itn>Newtmax-1) // if time step refused (Newton not converged)
 abort
 h = h0;
 t = t-dt;
 dt = dt/2
 restart = 1
 n = n -1;
else 
 restart = 0;  
end  
 

if (restart==0) // if time step accepted 


  h0 = h
     
 time(n) = time(n-1) + dt;
 timestep(n) = dt;
 

 dt = 1.2*dt;
 if (dt>dtmax)
  dt = dtmax;
 end
 if (t+dt>tf)
  dt = tf-t;
 end


// plot h and b 


////////////////////
// to be completed 
//////////////////// 
 
end // end if restart = 0

// plot Newton convergence
scf(1);
epsm = 1E-14;
plot(xn,log10(yn+epsm),'r-*')
clear xn;
clear yn;


end // end time loop ----------------------------------------------------------------



// computation of hs taking out erosions  and plot


////////////////////
//to be completed 
////////////////////

timer()
  
  




