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(So)
  s = Mo(So)/(Mo(So) + Mw(1-So));
endfunction  
  
  
  
function s = Fop(So)
  s = (Mop(So)*Mw(1-So) + Mo(So)*Mwp(1-So))/(Mo(So) + Mw(1-So))**2;
endfunction 



//data
Lx = 100;
Ly = 100;

//oil viscosity
muo=5E-3;

//water viscosity 
muw = 1.E-3

//residual oil saturation 
Sor = 0.3;

//irreducible water saturation
Swi = 0.1;

//productor well pressure 
pprod = 9.E+5;

//injector well pressure and injection concentration 
pinj = 10.E+5
Sinj = 0;

//well radius
rw = 0.1

//permeability (1 darcy) 
K = 1.E-12;
//K1 = 1.E-12;
//K2 = K1/100;

//porosity
phi = 0.1;




Day = 3600*24
tf = Day*360

Ndtmax = 10000



//space discretization: Nx and Ny (odd to have the well at the center)
Nx = 21;
Ny = Nx;
N = Nx*Ny;
dx = Lx/Nx;
dy = Ly/Ny;
Nint = Ny*(Nx-1) + Nx*(Ny-1);
Nbound = 2*(Nx+Ny);


//uniforme mesh of size dx dy 
// cells m = 1:N
for i=1:Nx
  Xi(i) = (i-1/2)*dx;
for j=1:Ny  
  Yj(j) = (j-1/2)*dy;
  
  m = i + (j-1)*Nx;
  
  volume(m) = dx*dy;
  X(m) = Xi(i);
  Y(m) = Yj(j);
end
end

//Interior faces: isig = 1:Nint
isig = 0;
//normal x 
for i=1:Nx-1
for j=1:Ny
  m1 = i+(j-1)*Nx;
  m2 = i+1+(j-1)*Nx;
  isig = isig + 1;
  cellint(isig,1) = m1;
  cellint(isig,2) = m2;
  surfaceint(isig) = dy;  
  Xint(isig) = (X(m1)+X(m2))/2;
  Yint(isig) = Y(m1);
end
end 
//normal y
for i=1:Nx
for j=1:Ny-1
  m1 = i+(j-1)*Nx;
  m2 = i+ j*Nx;
  isig = isig + 1;
  cellint(isig,1) = m1;
  cellint(isig,2) = m2;
  surfaceint(isig) = dx;
  Xint(isig) = X(m1);
  Yint(isig) = (Y(m1)+Y(m2))/2;
end
end 


//boundary faces x = 0
for j=1:Ny
  isig = j;
  i=1;
  m1 = i+(j-1)*Nx;
  cellbound(isig) = m1;
  surfacebound(isig) = dy;
  Xbound(isig) = 0;
  Ybound(isig) = Y(m1);
end

//boundary faces x = Lx
for j=1:Ny
  isig = j+Ny;
  i=Nx;
  m1 = i+(j-1)*Nx;
  cellbound(isig) = m1;
  surfacebound(isig) = dy;
  Xbound(isig) = Lx;
  Ybound(isig) = Y(m1);
end

//boundary faces y=0
for i=1:Nx
  isig = i+2*Ny;
  j=1;
  m1 = i+(j-1)*Nx;
  cellbound(isig) = m1;
  surfacebound(isig) = dx;
  Xbound(isig) = X(m1);
  Ybound(isig) = 0;
end

//boundary faces y=Ly
for i=1:Nx
  isig = i+Nx+2*Ny;
  j=Ny;
  m1 = i+(j-1)*Nx;
  cellbound(isig) = m1;
  surfacebound(isig) = dx;
  Xbound(isig) = X(m1);
  Ybound(isig) = Ly;
end


for m=1:N
  //permeabilities
  
  perm(m) = K;
  
  //  s = sin(20*X(m)/Lx)*sin(20*Y(m)/Ly)
   // if (s>0)
   //   perm(m) = K1;
   // else 
   //   perm(m) = K2;
   // end

end

//transmissibilities of interior faces 
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
  
  d1 = (Xint(isig)-X(m1))**2 + (Yint(isig)-Y(m1))**2;
  d1 = sqrt(d1);
  d2 = (Xint(isig)-X(m2))**2 + (Yint(isig)-Y(m2))**2;
  d2 = sqrt(d2);
  
  s1 = perm(m1)*surfaceint(isig)/d1;
  s2 = perm(m2)*surfaceint(isig)/d2;
  
  Tint(isig) = 1/(1/s1+1/s2);
end



//transmissibilities of boundary faces
for isig=1:Nbound
  m1 = cellbound(isig);
  d = (Xbound(isig)-X(m1))**2 + (Ybound(isig)-Y(m1))**2;
  d = sqrt(d);
  Tbound(isig) = perm(m1)*surfacebound(isig)/d;
end



//wells 

Nbwell = 5; //number of wells

r0 = 0.14*sqrt(dx**2+dy**2)
pi = 4*atan(1);

//1 productor well: at the center 

iw = int((Nx+1)/2) 
jw = int((Ny+1)/2)
mw = iw + (jw-1)*Nx
cellwell(1) = mw;
pwell(1) = pprod

// 4 injector wells: at the four corners

iw = 1
jw = 1
mw = iw + (jw-1)*Nx
cellwell(2) = mw;
pwell(2) = pinj

iw = 1
jw = Ny
mw = iw + (jw-1)*Nx
cellwell(3) = mw;
pwell(3) = pinj

iw = Nx
jw = Ny
mw = iw + (jw-1)*Nx
cellwell(4) = mw;
pwell(4) = pinj

iw = Nx
jw = 1
mw = iw + (jw-1)*Nx
cellwell(5) = mw;
pwell(5) = pinj


for iw=1:Nbwell
    mw = cellwell(iw)
    Twell(iw) = perm(mw)*2*pi/log(r0/rw)
end


//maximum of the fractional flow derivative: maxFop 
maxFop = 0;
for i=1:1000
  xx(i) = (i-1)/999;
  yy(i) = Fop(xx(i));
  if (yy(i)>maxFop)
    maxFop = yy(i);
  end
end



// initialization of the tracer 

t = 0;
n = 0;

//initial saturation S0 
for m=1:N
   S0(m) = 1-Swi;
end

//initial total mobilities MTint and MTbound 

// inner faces MTint
for isig=1:Nint
MTint(isig) = Mo(1-Swi) + Mw(Swi);
end


MTwell(1) = Mo(1-Swi) + Mw(Swi);
for iw=2:Nbwell
MTwell(iw) = Mw(1-Sinj) + Mo(Sinj);    
end



S = S0; 




//time loop ----------------------
while ( (t<tf)&(n<Ndtmax) )
// -------------------------------


n = n + 1


//----------------------------------

if ((int(n/10)*10==n)|(n==1)) then
//if (n==n) then

  
// pressure computation 
//Jacobian matrix and right hand side 
A=sparse([],[],[N,N]);
 
for m=1:N
  B(m,1) = 0;
end

 //Loop on interior faces
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
  
  s = Tint(isig)*MTint(isig);
  
  A(m1,m1) = A(m1,m1) + s;
  A(m2,m2) = A(m2,m2) + s; 
  A(m1,m2) = A(m1,m2) - s;
  A(m2,m1) = A(m2,m1) - s;
end



 //productor well: imposed pressure
 for iw = 1:Nbwell
 
 mw = cellwell(iw);
 s = Twell(iw)*MTwell(iw);
 A(mw,mw) = A(mw,mw) + s;
 B(mw,1) = B(mw,1) + s*pwell(iw); 
 
 end

 

 //solution
 p=A\B;
 


  
for i=1:Nx
for j=1:Ny  
  m = i + (j-1)*Nx; 
  pxy(i,j) = p(m)-pprod;
end  
end




scf(1);
xtitle('pressure ')
plot3d(Xi/100,Yj/100,pxy,40,80)
  

// Computation of the Velocities 
//Inner faces velocities 
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
  Vint(isig) = MTint(isig)*Tint(isig)*(p(m1)-p(m2));
end
//Well velocities 
for iw=1:Nbwell
mw = cellwell(iw);
Vwell(iw) = MTwell(iw)*Twell(iw)*(p(mw)-pwell(iw));
end

// Upwind cell at inner faces

for i=1:Nint
  m1 = cellint(i,1)
  m2 = cellint(i,2)
  if (Vint(i)<0)
      cellupwind(i) = m2
  else 
      cellupwind(i) = m1      
  end 
end


///////////////////////////////////////////////
// Max of the sum of the input velocities at each cell 
///////////////////////////////////////////////
for m=1:N
  Vinput(m)=0;
end
for isig=1:Nint
  m1 =  cellint(isig,1);
  m2 = cellint(isig,2);
  if (Vint(isig)<0)
    Vinput(m1) = Vinput(m1)-Vint(isig);
  else 
    Vinput(m2) = Vinput(m2)+Vint(isig);  
  end
end
for iw=1:Nbwell
  mw = cellwell(iw);
  if (Vwell(iw)<0) then
     Vinput(mw) = Vinput(mw) - Vwell(iw);   
  end
end

// computation of the maximum time step 

dtmax = min(phi*volume./Vinput)/maxFop

print(%io(2), 'dtmax in Day ',dtmax/Day)

// time step with CFL 1 
dt = dtmax
//dt = dtmax*3 // test CFL = 3 !  > unstable 

end

t = t + dt;

// saturation 

// loop on inner faces 

for i=1:Nint
  mup = cellupwind(i)
  m1 = cellint(i,1)
  m2 = cellint(i,2)
  
  flux = Fo(S0(mup))*Vint(i);
 
  S(m1) = S(m1) - flux*dt/phi/volume(m1);
  S(m2) = S(m2) + flux*dt/phi/volume(m2);
  
end

// loop on wells 
for iw=1:Nbwell
    mw = cellwell(iw)
    if (Vwell(iw)>=0)
      flux = Vwell(iw)*Fo(S0(mw))            
    else
      flux = Vwell(iw)*Fo(Sinj) 
    end
    
    S(mw) = S(mw) - flux*dt/phi/volume(mw);

end


//initial total mobilities MTint and MTbound 

// inner faces MTint
for isig=1:Nint
    
// upwind approximation 
//    
//    if (Vint(isig)>=0) 
//        MTint(isig) = Mo(S(m1)) + Mw(1-S(m1));
//    else 
//        MTint(isig) = Mo(S(m2)) + Mw(1-S(m2));
//    end
//
//
// centred approximation 

    MTint(isig) = Mo( (S(m1) + S(m2))/2 ) + Mw(1- (S(m1) + S(m2))/2 )

end

mw = cellwell(1)
MTwell(1) = Mo(S(mw)) + Mw(1-S(mw))
for iw=2:Nbwell
MTwell(iw) = Mw(1-Sinj) + Mo(Sinj);    
end


time(n) = t/Day // time 
//tracer  flow rate at the producer well 
mw = cellwell(1) 
flowrateprodoil(n) =  Vwell(1)*Fo(S0(mw))   


flowrateprodwater(n) =  Vwell(1)*(1-Fo(S0(mw)))   



// update of S0
S0 = S;


// plots every 10 time steps 
//2D plots
if ((int(n/10)*10==n)|(n==1)) then
for i=1:Nx
for j=1:Ny  
  m = i + (j-1)*Nx; 
  Sxy(i,j) = 1-S(m);
end  
end
scf(2);
xtitle('water saturation ')
plot3d(Xi/100,Yj/100,Sxy,40,80)
//diagonal plot
for i=1:Nx
    m = i + (i-1)*Nx
    Swdiag(i) = 1-S(m);
    Xdiag(i) = sqrt(X(m)**2+Y(m)**2)
end
scf(4)
xtitle(' diag plot of Sw ')
plot(Xdiag,Swdiag,'-r')
end
 
//////////////////////////////
end
//end of the time loop
//////////////////////////////

scf(3)
plot(time,flowrateprodoil,'-xr')
plot(time,flowrateprodwater,'-xb') 
  
  
scf(5)
xtitle(' diag plot of Sw ')
plot(Xdiag,Swdiag,'-xb')
  
