clear

//data
Lx = 100;
Ly = 100;

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

//injector well pressure and injection concentration 
pinj = 10.E+5
cinj = 1;

//well radius
rw = 0.1

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

//viscosity 
visco = 1.E-3

//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 = 11;
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;
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 

mw = int((Nx+1)/2) + (int((Ny+1)/2)-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




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


  
// 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);
  
  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);
 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) = 
end
//Well velocities 
for iw=1:Nbwell
mw = cellwell(iw);
//Vwell(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) = 
  else 
     // cellupwind(i) = 
  end 
end


// choice of the time step 

// initialization of the tracer 

t = 0;
n = 0;
c = zeros(N,1);
c0 = c;


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

t = t + dt;
n = n + 1

// computation of c at new time step 

// loop on inner faces 

for i=1:Nint

  
end

// loop on wells 
for iw=1:Nbwell


end
///////

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


// sum of the injector flow rates 
flowrateinj(n) = 0
for iw=2:Nbwell
  mw = cellwell(iw) 
  flowrateinj(n) =  flowrateinj(n) - Vwell(iw)*cinj     
end


// update of c0
c0 = c;


for i=1:Nx
for j=1:Ny  
  m = i + (j-1)*Nx; 
  cxy(i,j) = c(m);
end  
end
scf(2);
xtitle('tracer ')
plot3d(Xi/100,Yj/100,cxy,40,80)



 
//////////////////////////////
end
//end of the time loop
//////////////////////////////

scf(3)
plot(time,flowrateprod,'-xr')
plot(time,flowrateinj,'-xb') 
  
  

  
