clear






function s = Mo(So)
  if (So<Sor)
   s = 0;
  else 
   s = (So-Sor)**2/(1-Sor)**2/muo
  end
endfunction  
  
  
function s = Mw(Sw)
 if (Sw<Swi) 
  s = 0;
 else 
  s = (Sw-Swi)**2/(1-Swi)**2/muw
 end
endfunction

function s = Mop(So)
  if (So<Sor)
   s = 0;
  else 
   s = 2*(So-Sor)/(1-Sor)**2/muo
  end  
endfunction  
  
function s = Mwp(Sw)
 if (Sw<Swi) 
  s = 0;
 else 
  s = 2*(Sw-Swi)/(1-Swi)**2/muw
 end
endfunction

function s = Fo(S)
  s = Mo(S)/(Mo(S) + Mw(1-S));
endfunction  
  
  
  
function s = Fop(S)
  s = (Mop(S)*Mw(1-S) + Mo(S)*Mwp(1-S))/(Mo(S) + Mw(1-S))**2;
endfunction 



////// Residuals incompressible immiscible Two Phase darcy flow 

function [fluxint,fluxbound,upwind,R] = residu(p,S,p0,S0,dt)

//Inner Fluxes -\int_sigma Kgrad(p).n at face i oriented from cell 1 to cell 2
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);
  fluxint(i) = Tint(i)*(p(m1)-p(m2));
end

//Boundary fluxes 
for i=1:Nbound
    m = cellbound(i);
    fluxbound(i) = Tbound(i)*(p(m)-pbound(i));
end


//Upwind cell of the inner face i 
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);  
  if (fluxint(i)>=0)
   upwind(i) = m1;
  else
   upwind(i) = m2;
  end
end

//
//
// Oil and Water residuals in each cell 
//
//
// Accumulation term 
for m=1:N
  R(m,1) =  rhoo*phi*volume(m)/dt*( S(m) - S0(m) );
  R(m,2) = -rhow*phi*volume(m)/dt*( S(m) - S0(m) );
end

// Inner Fluxes
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);
  j = upwind(i);
  s1 = fluxint(i)*rhoo*Mo(S(j));
  s2 = fluxint(i)*rhow*Mw(1-S(j));
  R(m1,1) = R(m1,1) + s1;
  R(m2,1) = R(m2,1) - s1;
  R(m1,2) = R(m1,2) + s2;
  R(m2,2) = R(m2,2) - s2;  
end

// Boundary fluxes 
for i=1:Nbound
    m = cellbound(i);
    if (fluxbound(i)<0)
     //input boundary  
     R(m,1) =  R(m,1) + fluxbound(i)*rhoo*Mo(Sbound(i));
     R(m,2) =  R(m,2) + fluxbound(i)*rhow*Mw(1-Sbound(i));
    else 
     //Ouput water and oil from the cell m
     R(m,1) =  R(m,1) + fluxbound(i)*rhoo*Mo(S(m));  
     R(m,2) =  R(m,2) + fluxbound(i)*rhow*Mw(1-S(m));
    end
end

endfunction


//
//
//    Jacobian Matrix  
//
//
function A = Jacobian(p,S,fluxint,fluxbound,upwind,dt)

//Jacobian 
  A=sparse([],[],[2*N,2*N]); // sparse matrix format 
// A = zeros(2*N,2*N); // full matrix format
// 
//    Storage order: A(ii,jj)      ii=1,...,2*N,  jj=1,...,2*N
//
//    Equations: oil and water --> ii=2*i-1 (oil),2*i (water), i=1,...,N
//
//    Unknowns:   p and S=So   --> jj=2*j-1 (p),2*j (S), j=1,...,N 
//
//
//
// Accumulation terms
for m=1:N

  // Oil equation in cell m
  A(2*m-1,2*m-1) = 0; // derivative wrt p(m)
  A(2*m-1,2*m) = rhoo*phi*volume(m)/dt; // derivative wrt S(m)
  
  // Water equation in cell m 
  A(2*m,2*m-1) = 0; // derivative wrt p(m)
  A(2*m,2*m)   = -rhow*phi*volume(m)/dt; // derivative wrt S(m)
end

//Fluxes 
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);
  

j = upwind(i)

//derivatives wrt p(m1) and p(m2)
//
// oil equations in cells m1 and m2 
ss = Tint(i)*rhoo*Mo(S(j));  
A(2*m1-1,2*m1-1) = A(2*m1-1,2*m1-1) + ss;
A(2*m1-1,2*m2-1) = A(2*m1-1,2*m2-1) - ss;  
A(2*m2-1,2*m1-1) = A(2*m2-1,2*m1-1) - ss;  
A(2*m2-1,2*m2-1) = A(2*m2-1,2*m2-1) + ss;

// water equations  in cells m1 and m2
ss = Tint(i)*rhow*Mw(1-S(j));
A(2*m1,2*m1-1) = A(2*m1,2*m1-1) + ss;
A(2*m1,2*m2-1) = A(2*m1,2*m2-1) - ss;  
A(2*m2,2*m1-1) = A(2*m2,2*m1-1) - ss;  
A(2*m2,2*m2-1) = A(2*m2,2*m2-1) + ss;  


//derivatives wrt the upwind saturation S(j)
//
// Oil equations in cells m1 and m2 
ss = fluxint(i)*rhoo*Mop(S(j));  
A(2*m1-1,2*j) = A(2*m1-1,2*j) + ss;
A(2*m2-1,2*j) = A(2*m2-1,2*j) - ss;



// Water equations in cells m1 and m2 
ss = fluxint(i)*rhow*Mwp(1-S(j))
A(2*m1,2*j) = A(2*m1,2*j) - ss
A(2*m2,2*j) = A(2*m2,2*j) + ss

end

//
//
//Boundary fluxes 
//
for i=1:Nbound
    m = cellbound(i);
    if (fluxbound(i)<0)
     //input boundary  
     //derivatives wrt p(m) for the oil and water equations in cell m 
      A(2*m-1,2*m-1) = A(2*m-1,2*m-1) + Tbound(i)*rhoo*Mo(Sbound(i));   
      A(2*m,2*m-1) = A(2*m,2*m-1) + Tbound(i)*rhow*Mw(1-Sbound(i));     
     
    else 
     // output boundary 
     //derivatives wrt p(m) for the oil and water equations in cell m        
      A(2*m-1,2*m-1) = A(2*m-1,2*m-1) + Tbound(i)*rhoo*Mo(S(m));   
      A(2*m,2*m-1) = A(2*m,2*m-1) + Tbound(i)*rhow*Mw(1-S(m));    

     //derivatives wrt S(m) for the oil and water equations in cell m 
      A(2*m-1,2*m) = A(2*m-1,2*m) + fluxbound(i)*rhoo*Mop(S(m));   
      A(2*m,2*m) = A(2*m,2*m)     - fluxbound(i)*rhow*Mwp(1-S(m));
    end
end


//End of Jacobian computation  
 

endfunction







//data set   //////////////////////////////////////////////////

L=1000;
tf=3600*24*360*20;

//water viscosity
muw=1E-3;

//oil viscosity
muo=5E-3;

//permeability
K=1E-12;

//water density 
rhow = 1;

//oil density 
rhoo = 1; 

//porosity 
phi = 0.1;

//presssure at x=0: injector  
pinj = 1.1E+7;

//pressure at x=L: productor
pw=9.5E+6;

//Oil residual saturation
Sor = 0.3;

//Water irreductible saturation 
Swi = 0.1;

//Space discretization 
N=100;
Nint = N-1;
Nbound = 2;
dx=L/N;

//Time discretization
Ndtmax = 1000;
dtinit = 3600*24*5;
dtmax = 3600*24*5;

//Maximum number of Newton iterations
Newtmax = 25;
//Newton stopping criteria
eps=1.0E-6;
dSobj = 0.1;



// Mesh data structure  ////////////////////////////////////////

//Uniform 1D mesh of step dx 
// cells m = 1:N
for m=1:N
  volume(m) = dx;
  X(m) = (m-1/2)*dx;
end  

//Inner faces isig = 1:Nint
for isig=1:Nint
  cellint(isig,1) = isig;
  cellint(isig,2) = isig+1;
  surfaceint(isig) = 1;  
  Xint(isig) = isig*dx;
end

//boundary faces data structure
  cellbound(1) = 1;
  surfacebound(1) = 1;
  Xbound(1) = 0;
  
  cellbound(2) = N;
  surfacebound(2) = 1;
  Xbound(2) = L;


for m=1:N
  //permeability
    perm(m) = K;
end

//Inner face transmissivities
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);

  s1 = surfaceint(isig)*perm(m1)/abs(Xint(isig)-X(m1));
  s2 = surfaceint(isig)*perm(m2)/abs(Xint(isig)-X(m1));  
  
  Tint(isig) = 1/(1/s1+1/s2);
end
//Boundary face transmissivities
isig = 1;
m = cellbound(isig);
Tbound(1) = surfacebound(isig)*perm(m)/abs(Xbound(isig)-X(m));
isig = 2;
m = cellbound(isig);
Tbound(2) = surfacebound(isig)*perm(m)/abs(Xbound(isig)-X(m));


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


//pressure aand input oil saturation at the boundaries 
pbound(1) = pinj;
pbound(2) = pw;
Sbound(1) = 0; // water injection 
Sbound(2) = 0; // not used, output boundary

// initial pressure p and oil saturation S 

for i=1:N
  p0(i) = (pbound(1)+pbound(2))/2;  //Initial pressure 
  S0(i)=1-Swi;    //Initial oil saturation
end


// initialtime step
dt = dtinit;
temps(1) = 0;
pasdetemps(1) = dt;
//initial oil and water flow rates
flowrateTo(1) = 0;
flowrateTw(1) = 0;
flowrateo(1) = 0;
flowratew(1) = 0;

//static pressure for plotting
for m=1:N
  pstatic(m) = X(m)/L*pw+(L-X(m))/L*pinj;
end


//////// TIME LOOP ////////////////////////////////////////////
t = 0;
n = 0;
newt = 1;
timer()
while ((t<tf)&(n<Ndtmax))
n = n + 1;
t = t + dt
////////////////////////////////////////////////////////////////



//initialization of Newton 
p = p0;
S = S0;

[fluxint,fluxbound,upwind,R] = residu(p,S,p0,S0,dt)


//Initial residual norm 
normR = norm(R(:,1)) + norm(R(:,2));
normR0 = normR;


/////// Newton LOOP  ////////////////////////////////
itn = 1; 
xn(itn) = newt;
yn(itn) = 1; // xn, yn for plot of Newton convergence 
// ----------------------------------------------
while ((normR/normR0>eps)&(itn<Newtmax)) do 
// ----------------------------------------------  
itn = itn + 1;
newt = newt + 1;
xn(itn) = newt;

//Jacobian 
A = Jacobian(p,S,fluxint,fluxbound,upwind,dt)

 
//Right hand side = -Residual 
for i=1:N
  B(2*i-1,1)   = -R(i,1);
  B(2*i,1)     = -R(i,2);
end

//solution dU 
dU = A\B;
//clear A;

//max dS
dSmax=0;
for i=1:N
dSmax = max(dSmax,abs(dU(2*i)));
end
alpha = min(1,dSobj/dSmax);
alpha = 1

//incrementation of Newton 
for i=1:N 
 p(i) = p(i) + dU(2*i-1)*alpha;  
 S(i) = S(i) + dU(2*i)*alpha;
end


//residual 
[fluxint,fluxbound,upwind,R] = residu(p,S,p0,S0,dt)



//norm of the residual  
normR = norm(R(:,1)) + norm(R(:,2));
yn(itn) = normR/normR0; // for plot of Newton convergence



end
//  End Newton LOOP /////////////////////////////////////////


////////// if Newton not converged --> restart of the time step  
if (itn>Newtmax-1) 
 S = S0;
 p = p0;
 t = t-dt;
 dt = dt/2
 restart = 1
 n = n -1;
else 
 restart = 0;  
end  
/////////////

 
//Time step management and oil and water flow rates
 if (restart==0) 
 temps(n+1) = temps(n) + dt;
 pasdetemps(n+1) = dt;
 if (fluxbound(2)>0)
  flowrateo(n+1) = fluxbound(2)*Mo(S(N));  
  flowratew(n+1) = fluxbound(2)*Mw(1-S(N));
  
  flowrateTo(n+1) = flowrateTo(n) + flowrateo(n+1)*dt;  
  flowrateTw(n+1) = flowrateTw(n) + flowratew(n+1)*dt;
 end
 S0=S;
 p0=p,
 dt = 1.2*dt;
 if (dt>dtmax)
  dt = dtmax;
 end
 if (t+dt>tf)
  dt = tf-t;
 end

//// plots ////////

 //Oil saturation and pressure plots 
 if (int(n/1)*1==n)
  scf(12);
  title(' water saturation ')
  plot (X,1-S)
  scf(13);
  title(' pressure ')
  plot (X,p-pstatic)  
 end
end

//Newton convergence plot at each time step
scf(15);
title(' Newton convergence logarithm of the relative resiudal')
xtitle('Total number of Newton iterations')
plot(xn,log10(yn),'r-*')
clear xn;
clear yn;


end
//End time LOOP //////////////////////////////////////////////////////////////////

timer()

//Oil and water flow rates at x=L
scf(7);
title(' Oil and water flow rates ')
xlabel('Time in years')
plot(temps/3600/24/360,flowrateo,'-r')
plot(temps/3600/24/360,flowratew,'-b')

scf(17);
title(' Oil and Water total flow rate ')
xlabel('Time in years')
plot(temps/3600/24/360,flowrateTo,'-r')
plot(temps/3600/24/360,flowrateTw,'-b')

//Time steps plot 
scf(18);
title(' Time step function of time ')
xlabel('Time in years')
plot(temps/3600/24/360,pasdetemps/3600/24/360,'-*')


