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 = 1000; //domain (0,Lx)
tf = 3600*24*30*12*10;   //final simulation time

//production bottom well pressure
pw = 9.5E+6;

//injection well pressure
pinj = 1.1E+7;

//oil viscosity
muo=5E-3;

//water viscosity 
muw = 1.E-3

//residual oil saturation 
Sor = 0.3;

//irreducible water saturation
Swi = 0.1;

//reservoir permeability
K=1E-12;

//reservoir porosity
phi = 0.1;


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

dx=Lx/N;


/////////////////////////
//MESH DATA STRUCTURE 
/////////////////////////

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

//inner faces data structure: 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) = Lx;


//////////////////////////////////////////////////
//INNER AND BOUNDARY FACES TRANSMISSIVITIES
//////////////////////////////////////////////////

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

//inner faces 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 faces transmissivities
for isig=1:Nbound
  m = cellbound(isig);
  Tbound(isig) = surfacebound(isig)*perm(m)/abs(Xbound(isig)-X(m));
end


//////
//pressure boundary conditions 
//////

pbound(1) = pinj;
pbound(2) = pw;



//time
time(1) = 0;
//flowrates 
oilflowrateT(1) = 0;
oilflowrate(1) = 0;
waterflowrateT(1) = 0;
waterflowrate(1) = 0;


//fractional flow Fo function of oil saturation
for i=1:1000
  xx(i) = (i-1)/999;
  yy(i) = Fo(xx(i));
end
scf(1);
xtitle('Fractional flow Fo and its derivative') ;
plot(xx,yy)

//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
plot(xx,yy)


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


///////////////////////////////////////////////////
//INITIALIZATION: saturation S0 and total mobilities MTint and MTbound
///////////////////////////////////////////////////


n = 1;
t = 0;

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

//initial total mobilities MTint and MTbound 

// inner faces MTint
for isig=1:Nint
//     MTint(isig) = ;
end

// left end total mobility (water injection)
// MTbound(1) = ;

// right end total mobility (producer)
// MTbound(2) = ;


S = S0;


//////////////////////////////////////
//TIME LOOP
while ( (t<tf)&(n<10000) )
//
//////////////////////////////////////



////////////////////////////////
///////////////////////////////
//PRESSURE COMPUTATION p
///////////////////////////////
///////////////////////////////

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

//inner faces loop
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
  
end

//boundary faces loop
for isig=1:Nbound
  m = cellbound(isig);
  
  
end

 //pressure solution 
p=A\B;


/////////////////////////////////////////////////// 
///////////////////////////////////////////////////
//TOTAL VELOCITIES COMPUTATION  VTint and VTbound
///////////////////////////////////////////////////
///////////////////////////////////////////////////

//total velocities at inner faces VTint 
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
//  VTint(isig) = ;
end

//total velocities at boundary faces VTbound
for isig=1:Nbound
  m = cellbound(isig);
//  VTbound(isig) = ;
end


////////////////////////////////////////
////////////////////////////////////////
// CFL COMPUTATION FOR THE TIME STEP dt
////////////////////////////////////////
////////////////////////////////////////


//time step

dt = 5*3600*24;

//plot of the time step function of time 
scf(2);
xtitle('Time step function of time') ;
plot(t/(3600*24),dt/(3600*24),'*')


////////////////////////////////////////////////
///////////////////////////////////////////////
//COMPUTATION OF THE SATURATION S
///////////////////////////////////////////////
////////////////////////////////////////////////




//inner faces loop
for isig=1:Nint
 
 
 
end

//boundary faces loop 
for isig=1:Nbound


end 


/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
//COMPUTATION OF TOTAL MOBILITIES MTint and MTbound  
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////

//inner faces total mobilities MTint 
for isig=1:Nint
    m1 = cellint(isig,1);
    m2 = cellint(isig,2);

end

//boundary faces total mobilities MTbound 
for isig=1:Nbound
    m = cellbound(isig);

end

 
////////////////// 
//oil and water flowrates and total flowrates at the producer
////////////////// 
 
 time(n+1) = time(n) + dt;
 m = cellbound(2);
 oilflowrate(n+1) =  VTbound(2)*Fo(S0(m));
 oilflowrateT(n+1) = oilflowrateT(n) + oilflowrate(n+1)*dt;
 
 waterflowrate(n+1) =  VTbound(2)*(1-Fo(S0(m)));
 waterflowrateT(n+1) = waterflowrateT(n) + waterflowrate(n+1)*dt;


//////////////////////////////////////
///// update of time and of saturation at previous time step
//////////////////////////////////////

 t = t+dt;
 n = n + 1;
 
 S0 = S;
 
///////////////////////////
/////////////////////////// 
 
 
//////////////////////////// 
//PLOTS AT EACH TIME STEP 
////////////////////////////

if ((n==1)|(int(n/10)*10==n))
  scf(3);
  xtitle('OIL SATURATION') ;
  plot(X,S,'r-')
  
  scf(4);
  xtitle('WATER SATURATION') ;
  plot(X,1-S)
  

  
  scf(5);
  xtitle('OVER PRESSURE') 
  plot(X,p-pstatic)
end

//////////////////////////////////
// END OF THE TIME STEP LOOP
end
//////////////////////////////////



scf(6);
xtitle('o and w total flowrates') ;
plot(time,oilflowrateT,'r-*')
plot(time,waterflowrateT, '-*')


scf(7);
xtitle('o and w flowrates') ;
plot(time,oilflowrate,'r-*')
plot(time,waterflowrate, '-*')



