clear



function s = f_cbar(u)
  p0 = 7E+6;
  p1 = 1.1E+7;
 // s = 0.02*(u-p1)/(p0-p1) + 0.06*(u-p0)/(p1-p0) ;
   s = 0.1*(u-p1)/(p0-p1) + 0.2*(u-p0)/(p1-p0) ;
endfunction

function s = dfsdp_cbar(u)
  p0 = 7E+6;
  p1 = 1.1E+7;
 // s = 0.02/(p0-p1) + 0.06/(p1-p0) ;  
   s = 0.1/(p0-p1) + 0.2/(p1-p0) ;  
endfunction


function s = Mo(So)
 if (So<Soi) 
  s = 0 ;
 else 
  s = (So-Soi)/(1-Soi)/muo ;
 end
endfunction

function s = Mg(Sg)
  if (Sg<Sgr)
   s = 0 ;
  else 
   s = (Sg-Sgr)/(1-Sgr)/mug ;
  end
endfunction  

function s = Mop(So)
  if (So<Soi)
   s = 0;
  else 
   s = 1/(1-Soi)/muo ;
  end  
endfunction  
  
function s = Mgp(Sg)
 if (Sg<Sgr) 
  s = 0;
 else 
  s = 1/(1-Sgr)/mug ;
 end
endfunction



function s = f_rhoo(u,c)
  p0 = 7E+6;
  p1 = 1.1E+7;
  s = 800*(u-p1)/(p0-p1) + 1000*(u-p0)/(p1-p0) + 800*c ;
endfunction

function s = dfsdp_rhoo(u,c)
  p0 = 7E+6;
  p1 = 1.1E+7;
  s = 800/(p0-p1) + 1000/(p1-p0) ;
endfunction

function s = dfsdc_rhoo(u,c)
  s = 800 ;
endfunction

function s = f_rhog(u)
  p0 = 7E+6;
  p1 = 1.1E+7;
  s = 300*(u-p1)/(p0-p1) + 700*(u-p0)/(p1-p0);
endfunction

function s = dfsdp_rhog(u)
  p0 = 7E+6;
  p1 = 1.1E+7;
  s = 300/(p0-p1) + 700/(p1-p0);  
endfunction





function [rhoo,rhog,dprhoo,dcrhoo,dprhog,dpcbar,fluxint,fluxwell,upwind,R] = residu(p,S,C,DIPH,p0,S0,C0,dt)


for m=1:N
  
  rhoo(m) = f_rhoo(p(m),C(m));
  rhog(m) = f_rhog(p(m));
  dpcbar(m) = dfsdp_cbar(p(m));
  dprhog(m) = dfsdp_rhog(p(m));
  dcrhoo(m) = dfsdc_rhoo(p(m),C(m))
  
  if (DIPH(m)==0)
     dprhoo(m) = dfsdp_rhoo(p(m),C(m));
  else 
     dprhoo(m) = dfsdp_rhoo(p(m),C(m)) + dfsdc_rhoo(p(m),C(m))*dpcbar(m);
  end

end


//flux phase alpha -K(grad (p + rhoalpha*g*Y) ) aux aretes i oriente de la cell 1 vers la cell 2
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);
  fluxint(i,1) = Tint(i)*(p(m1)-p(m2) + grav*(rhoo0(m1)+rhoo0(m2))/2*(Y(m1)- Y(m2)));
  fluxint(i,2) = Tint(i)*(p(m1)-p(m2) + grav*(rhog0(m1)+rhog0(m2))/2*(Y(m1)- Y(m2))); 
end

//flux well
for i=1:Nwell
 fluxwell(i) = Twell(i)*( p(cellwell(i)) - pressionwell(i) );
end


//cell upwind pour la phase alpha et pour l'arete i interieure 
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);  
  
  if (fluxint(i,1)>=0)
   upwind(i,1) = m1;
  else
   upwind(i,1) = m2;
  end

  if (fluxint(i,2)>=0)
   upwind(i,2) = m1;
  else
   upwind(i,2) = m2;
  end

end

//calcul des residus 

//boucle sur les cells: accumulation Lourd et Leger
for m=1:N
  R(m,1) =  phi*volume(m)/dt*( rhoo(m)*(1-C(m))*(1-S(m)) - rhoo0(m)*(1-C0(m))*(1-S0(m)) ) ;
  R(m,2) =  phi*volume(m)/dt*( rhoo(m)*C(m)*(1-S(m)) - rhoo0(m)*C0(m)*(1-S0(m))  );
  R(m,2) =  R(m,2) + phi*volume(m)/dt*( rhog(m)*S(m) - rhog0(m)*S0(m) );
end

//boucle sur les aretes interieures 
for i=1:Nint
  m1 = cellint(i,1);
  m2 = cellint(i,2);  
  mo = upwind(i,1);
  mg = upwind(i,2); 
  s1 = fluxint(i,1)*rhoo(mo)*Mo(1-S(mo))*(1-C(mo));
  s2 = fluxint(i,1)*rhoo(mo)*Mo(1-S(mo))*C(mo) + fluxint(i,2)*rhog(mg)*Mg(S(mg)) ;
  
  R(m1,1) = R(m1,1) + s1;
  R(m1,2) = R(m1,2) + s2;
  R(m2,1) = R(m2,1) - s1;
  R(m2,2) = R(m2,2) - s2;  
end



//well producteur 
m = cellwell(1) ;
if (fluxwell(1)>=0)
   R(m,1) =  R(m,1) + fluxwell(1)*rhoo(m)*Mo(1-S(m))*(1-C(m));  
   R(m,2) =  R(m,2) + fluxwell(1)*rhoo(m)*Mo(1-S(m))*C(m);
   R(m,2) = R(m,2) + fluxwell(1)*rhog(m)*Mg(S(m));
end 


//well injecteur 
for i=2:Nwell
 m = cellwell(i);
 if (fluxwell(i)<0)
   R(m,2) =  R(m,2) + fluxwell(i)*f_rhog(pressionwell(i))*Mg(1);
 end
end

endfunction



function A = jacobienne(upwind,p,S,C,DIPH,dt,rhoo,rhog,dprhoo,dcrhoo,dprhog,dpcbar,fluxint,fluxwell)
    
    
//calcul de la Jacobienne
A=sparse([],[],[2*N,2*N]);

// Ordre: 
//  - équation composant lourd (impair) puis composant leger (pair)
//  - Inconnue p (pair) et C=C_leger dans l'huile en monophasique et Sg en diphasique (impair) 

//boucle sur les cells: terme d'accumulation
for i=1:N
  //Lourd - p
  A(2*i-1,2*i-1) = phi*volume(i)/dt*dprhoo(i)*(1-C(i))*(1-S(i))- DIPH(i)*phi*volume(i)/dt*rhoo(i)*dpcbar(i)*(1-S(i));
  //Lourd - C ou S
  A(2*i-1,2*i) = (phi*volume(i)/dt)*( (1-DIPH(i))*dcrhoo(i)*(1-S(i))*(1-C(i)) - (1-DIPH(i))*(1-S(i))*rhoo(i) - DIPH(i)*(1-C(i))*rhoo(i) ) ;
  //Leger - p
  A(2*i,2*i-1) = phi*volume(i)/dt*dprhog(i)*S(i) + phi*volume(i)/dt*dprhoo(i)*C(i)*(1-S(i)) + DIPH(i)*phi*volume(i)/dt*rhoo(i)*dpcbar(i)*(1-S(i));
  //Leger - C ou S
  A(2*i,2*i)   = (phi*volume(i)/dt)*( (1-DIPH(i))*dcrhoo(i)*C(i)*(1-S(i)) + (1-DIPH(i))*rhoo(i)*(1-S(i)) + DIPH(i)*(-rhoo(i)*C(i)+rhog(i)) );
end

//boucle sur les aretes interieures 
for i=1:Nint
  
  m1 = cellint(i,1);
  m2 = cellint(i,2);
  mo = upwind(i,1);
  mg = upwind(i,2);

//derivees par rapport aux pressions p(m1), p(m2) dans les potentiels

ss = Tint(i)*rhoo(mo)*Mo(1-S(mo))*(1-C(mo));  
//Lourd m1 fluxint p(m1)
A(2*m1-1,2*m1-1) = A(2*m1-1,2*m1-1) + ss;
//Lourd m1 fluxint p(m2)
A(2*m1-1,2*m2-1) = A(2*m1-1,2*m2-1) - ss;  
//Lourd m2 fluxint p(m1)
A(2*m2-1,2*m1-1) = A(2*m2-1,2*m1-1) - ss;  
//Lourd m2 fluxint p(m2)
A(2*m2-1,2*m2-1) = A(2*m2-1,2*m2-1) + ss;

ss = Tint(i)*rhoo(mo)*Mo(1-S(mo))*C(mo) + Tint(i)*rhog(mg)*Mg(S(mg)) ;
//Leger m1 fluxint p(m1)
A(2*m1,2*m1-1) = A(2*m1,2*m1-1) + ss;
//Leger m1 fluxint p(m2)
A(2*m1,2*m2-1) = A(2*m1,2*m2-1) - ss;  
//Leger m2 fluxint p(m1)
A(2*m2,2*m1-1) = A(2*m2,2*m1-1) - ss;  
//Leger m2 fluxint p(m2)
A(2*m2,2*m2-1) = A(2*m2,2*m2-1) + ss;  

//derivees par rapport a p(mo) dans rhoo eq Lourd
ss = fluxint(i,1)*dprhoo(mo)*Mo(1-S(mo))*(1-C(mo));  
//Lourd m1
A(2*m1-1,2*mo-1) = A(2*m1-1,2*mo-1) + ss;
//Lourd m2 
A(2*m2-1,2*mo-1) = A(2*m2-1,2*mo-1) - ss;


//derivees par rapport a C(mo) dans rhoo eq Lourd MONO
ss = (1-DIPH(mo))*fluxint(i,1)*dcrhoo(mo)*Mo(1-S(mo))*(1-C(mo));  
//Lourd m1
A(2*m1-1,2*mo) = A(2*m1-1,2*mo) + ss;
//Lourd m2 
A(2*m2-1,2*mo) = A(2*m2-1,2*mo) - ss;
//derivees par rapport a p(mo) dans (1-C) eq Lourd en DIPH=1
ss = -DIPH(mo)*fluxint(i,1)*rhoo(mo)*Mo(1-S(mo))*dpcbar(mo);  
//Lourd m1
A(2*m1-1,2*mo-1) = A(2*m1-1,2*mo-1) + ss;
//Lourd m2 
A(2*m2-1,2*mo-1) = A(2*m2-1,2*mo-1) - ss;
//derivees par rapport a p(mo) dans rhoo eq Leger
ss = fluxint(i,1)*dprhoo(mo)*Mo(1-S(mo))*C(mo);
//Leger m1 - phase huile
A(2*m1,2*mo-1) = A(2*m1,2*mo-1) + ss;
//Leger m2  - phase huile
A(2*m2,2*mo-1) = A(2*m2,2*mo-1) - ss;

//derivees par rapport a C(mo) dans rhoo eq Leger MONO
ss = (1-DIPH(mo))*fluxint(i,1)*dcrhoo(mo)*Mo(1-S(mo))*C(mo);
//Leger m1 - phase huile
A(2*m1,2*mo) = A(2*m1,2*mo) + ss;
//Leger m2  - phase huile
A(2*m2,2*mo) = A(2*m2,2*mo) - ss;


//derivees par rapport a p(mo) dans C eq Leger en DIPH=1
ss = DIPH(mo)*fluxint(i,1)*rhoo(mo)*Mo(1-S(mo))*dpcbar(mo);
//Leger m1 - phase huile
A(2*m1,2*mo-1) = A(2*m1,2*mo-1) + ss;
//Leger m2  - phase huile
A(2*m2,2*mo-1) = A(2*m2,2*mo-1) - ss;

//derivees par rapport a p(mg) dans rhog eq Leger
ss = fluxint(i,2)*dprhog(mg)*Mg(S(mg)) ;
//Leger m1 - phase gaz
A(2*m1,2*mg-1) = A(2*m1,2*mg-1) + ss;
//Leger m2 - phase gaz
A(2*m2,2*mg-1) = A(2*m2,2*mg-1) - ss ;


//derivees eq Lourd dans (1-C)*Mo par rapport a C en mono, S en diph
ss = fluxint(i,1)*rhoo(mo)*( - (1-DIPH(mo))*Mo(1-S(mo)) - DIPH(mo)*(1-C(mo))*Mop(1-S(mo)) );  
//Lourd m1
A(2*m1-1,2*mo) = A(2*m1-1,2*mo) + ss;
//Lourd m2 
A(2*m2-1,2*mo) = A(2*m2-1,2*mo) - ss;

//derivees eq Leger terme C*Mo  par rapport a C en mono, S en diph
ss = fluxint(i,1)*rhoo(mo)*((1-DIPH(mo))*Mo(1-S(mo)) - DIPH(mo)*C(mo)*Mop(1-S(mo))) ;
//Leger m1 - phase huile
A(2*m1,2*mo) = A(2*m1,2*mo) + ss;
//Leger m2  - phase huile
A(2*m2,2*mo) = A(2*m2,2*mo) - ss;

//derivees eq Leger terme Mg par rapport a C en mono, S en diph
ss = DIPH(mg)*fluxint(i,2)*rhog(mg)*Mgp(S(mg)) ;
//Leger m1 - phase gaz
A(2*m1,2*mg) = A(2*m1,2*mg) + ss;
//Leger m2 - phase gaz
A(2*m2,2*mg) = A(2*m2,2*mg) - ss ;

end




//well producteur 
m1 = cellwell(1) ;
if (fluxwell(1)>=0)
  
  //derivees par rapport a p(m1) dans les potentiels 
  A(2*m1-1,2*m1-1) = A(2*m1-1,2*m1-1) + Twell(1)*rhoo(m1)*Mo(1-S(m1))*(1-C(m1));  
  A(2*m1,2*m1-1)   = A(2*m1,2*m1-1)   + Twell(1)*( rhoo(m1)*Mo(1-S(m1))*C(m1) + rhog(m1)*Mg(S(m1)) ) ;  
 
  
//derivees par rapport a p(m1) dans rhoo et rhog
  A(2*m1-1,2*m1-1) = A(2*m1-1,2*m1-1) + fluxwell(1)*dprhoo(m1)*Mo(1-S(m1))*(1-C(m1)) ;   
  A(2*m1,2*m1-1)   = A(2*m1,2*m1-1)   + fluxwell(1)*dprhoo(m1)*Mo(1-S(m1))*C(m1) ; 
  A(2*m1,2*m1-1) = A(2*m1,2*m1-1)     + fluxwell(1)*dprhog(m1)*Mg(S(m1)) ;     
  
  
  //derivees par rapport a C(m1) dans rhoo en MONO
  A(2*m1-1,2*m1) = A(2*m1-1,2*m1) + (1-DIPH(m1))*fluxwell(1)*dcrhoo(m1)*Mo(1-S(m1))*(1-C(m1)) ;   
  A(2*m1,2*m1)   = A(2*m1,2*m1)   + (1-DIPH(m1))*fluxwell(1)*dcrhoo(m1)*Mo(1-S(m1))*C(m1) ; 
  
//derivees par rapport a p(m1) dans C et 1-C en diphasique
  A(2*m1-1,2*m1-1) = A(2*m1-1,2*m1-1) - DIPH(m1)*fluxwell(1)*rhoo(m1)*Mo(1-S(m1))*dpcbar(m1) ;   
  A(2*m1,2*m1-1)   = A(2*m1,2*m1-1)   + DIPH(m1)*fluxwell(1)*rhoo(m1)*Mo(1-S(m1))*dpcbar(m1) ; 
    
  
  
  //derivees par rapport a C(m1) en MONO et S(m1) en DIPH
  A(2*m1-1,2*m1) = A(2*m1-1,2*m1) + fluxwell(1)*rhoo(m1)*( -(1-DIPH(m1))*Mo(1-S(m1)) - DIPH(m1)*(1-C(m1))*Mop(1-S(m1)) );   
  A(2*m1,2*m1)   = A(2*m1,2*m1)   + fluxwell(1)*rhoo(m1)*( (1-DIPH(m1))*Mo(1-S(m1)) - DIPH(m1)*C(m1)*Mop(1-S(m1)) ); 
  A(2*m1,2*m1) = A(2*m1,2*m1)     + DIPH(m1)*fluxwell(1)*rhog(m1)*Mgp(S(m1));     
  
end 

//well injecteur 
for i=2:Nwell
 m1 = cellwell(i);
 if (fluxwell(i)<0)
  //Leger entrant dans le reservoir
  //derivees par rapport a p(m1) 
 A(2*m1,2*m1-1)   = A(2*m1,2*m1-1)   + Twell(i)*f_rhog(pressionwell(i))*Mg(1);  
 end
end

//fin calcul jacobienne 

endfunction




//-----------------------------------------------------//
// Cas 1D miscible, diphasique (O-G) BO avec composant leger dans le gaz et l'huile 
//  et composant lourd dans l'huile seule
// un well prod au centre, 2 well inj gaz aux deux bords      //
//
// C = concentration du leger dans l'huile 
//
//  rhoo(p,C)  rhog(p)
//
//  muo et mug constantes 
//
//  Cbar(p): concentration a l'equilibre 
//                                                     //
//-----------------------------------------------------//

//--------------------------------------------------
//                Données d'entrées
//--------------------------------------------------

// taille réservoir 
Lx= 100. ;
Ly = 100. ;
Lz = 1 ;
Nx = 201 ;
Ny = 1 ;

// temps simulation
tf = 3600*24*50 ;

//viscosite 
muo = 3E-4 ;
mug = 3E-5 ;


//points limites saturations
Soi = 0.3 ;
Sgr = 0.1 ;

//permeabilite
K=1E-14 ;

//porosite 
phi = 0.1 ;

//gravite
//grav = 9.81 ;
grav = 0 ;

//pression inj ;
pinj = 1.1E+7 ;

//pression prod
pprod = 7E+6 ;

//pression initiale
pinit = 1E+7 ;

//rayon well
rw = 0.1 ;

//--------------------------------------------------
//  Newton : max itérations + critère arret
//--------------------------------------------------

Newtmax = 25 ; 
eps=1.0E-10 ;


//--------------------------------------------------
//  Discrétisation en temps
//--------------------------------------------------

Ndtmax = 10000 ;
dtinit = 3600*24/10 ;
dtmax = 3600*24 ;


//--------------------------------------------------
//     Structure de donnees du maillage :
//
// discrétisation, info sur maillage, transmissivité
//--------------------------------------------------

//discretisation en espace
N = Nx*Ny ;
dx = Lx/Nx ;
dy = Ly/Ny ;
Nint = Ny*(Nx-1) + Nx*(Ny-1) ;
Nbord = 2*Ny ;

//maillage uniforme 2D de pas 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*Lz ;
  X(m) = Xi(i) ;
  Y(m) = Yj(j) ;
  
end
end

//structure de donnees aretes interieures: isig = 1:Nint
isig = 0 ;
//aretes verticales 
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 ;
  surface(isig) = dy*Lz ;  
  Xsig(isig) = (X(m1)+X(m2))/2 ;
  Ysig(isig) = Y(m1) ;
end
end 

//aretes horizontales
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 ;
  surface(isig) = dx*Lz ;
  Xsig(isig) = X(m1) ;
  Ysig(isig) = (Y(m1)+Y(m2))/2 ;
end
end 

for m=1:N
  //permeabilites
    if ((X(m)>6*Lx/10)&(X(m)<8*Lx/10))
     //perm(m) = K/100 ;
     perm(m) = K ;
    else
     perm(m) = K ;    
    end
end

//transmissibilites des aretes interieures 
for isig=1:Nint
  m1 = cellint(isig,1) ;
  m2 = cellint(isig,2) ;
  
  d1 = (Xsig(isig)-X(m1))**2 + (Ysig(isig)-Y(m1))**2 ;
  d1 = sqrt(d1) ;
  d2 = (Xsig(isig)-X(m2))**2 + (Ysig(isig)-Y(m2))**2 ;
  d2 = sqrt(d2) ;
  
  s1 = perm(m1)*surface(isig)/d1 ;
  s2 = perm(m2)*surface(isig)/d2 ;
  
  Tint(isig) = 1/(1/s1+1/s2) ;
end

//transmissibilite des well vertical

//1 well producteur et 2 well injecteurs
Nwell = 3;

//indice de productivite des well
pi = 4*atan(1);

//well producteur 
ifond = int((Nx+1)/2);
jfond = int((Ny+1)/2);
cellwell(1) = ifond+(jfond-1)*Nx ;
r0 = 0.28*sqrt(dx**2+dy**2)/2;
Twell(1) = perm(cellwell(1))*Lz*2*pi/log(r0/rw);
pressionwell(1) = pprod;

//well injecteurs 
ifond = 1 ;
jfond = 1 ;
cellwell(2) = ifond+(jfond-1)*Nx ;
r0 = 0.28*sqrt(dx**2+dy**2)/2;
Twell(2) = perm(cellwell(2))*Lz*2*pi/log(r0/rw);
pressionwell(2) = pinj;

ifond = Nx ;
jfond = 1 ;
cellwell(3) = ifond+(jfond-1)*Nx ;
r0 = 0.28*sqrt(dx**2+dy**2)/2;
Twell(3) = perm(cellwell(3))*Lz*2*pi/log(r0/rw);
pressionwell(3) = pinj;




//--------------------------------------------------
//                 Initialisation
//--------------------------------------------------

// Init des 3 inconnues : 
//       p : pression
//       S : saturation en gaz
//       C : concentration en composant leger dans la phase huile

for m=1:N
   S0(m) = 0 ;
   p0(m) = pinit ;
   C0(m) = f_cbar(p0(m));
   // Constante pour nombre de phase
   //DIPH = 0 si monophasique huile et DIPH = 1 si diphasique huile-gaz
   DIPH_0(m) = 0 ;   
   rhoo0(m) = f_rhoo(p0(m),C0(m));
   rhog0(m) = f_rhog(p0(m));
end

//temps
dt = dtinit ;
temps(1) = 0 ;
pasdetemps(1) = dt ;

//debits et cumul au well prod
debitg(1) = 0;
debito(1) = 0;
debitTg(1) = 0;
debitTo(1) = 0;

//--------------------------------------------------
//                   Itération
//--------------------------------------------------

//boucle en temps -----------------------
t = 0 ;
n = 1 ;
newt = 1 ;
timer()
while ((t<tf*0.999999999)&(n<Ndtmax))
n = n + 1 ;
t = t + dt ;
// boucle en temps ----------------------

//init Newton 
p = p0 ;
C = C0 ;
S = S0 ;
DIPH = DIPH_0 ;

[rhoo,rhog,dprhoo,dcrhoo,dprhog,dpcbar,fluxint,fluxwell,upwind,R] = residu(p,S,C,DIPH,p0,S0,C0,dt)

//calcul de la norme du residu 
normR = 0 ;
for i=1:N
  normR = normR + R(i,1)**2 + R(i,2)**2 ;
end
normR = sqrt(normR) ;
normR0 = normR ;

//boucle Newton 
itn = 1 ;
xn(itn) = newt ;
yn(itn) = 1 ;
alphan(itn) = 1;
for i=1:N
    idpl(i) = 0; // indicateur de passage au point limite Sg = 1-Sor 
end
// ----------------------------------------------
while ((normR/normR0>eps)&(itn<Newtmax)&(normR>1.0E-10)) do 
// ----------------------------------------------  
itn = itn + 1;
newt = newt + 1;
xn(itn) = newt;

A = jacobienne(upwind,p,S,C,DIPH,dt,rhoo,rhog,dprhoo,dcrhoo,dprhog,dpcbar,fluxint,fluxwell)
 
//second membre
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;


for i=1:N
  dSC(i) = dU(2*i);
end
//scf(20)
//title('dSC')
//plot(Xi,dSC)

// global relaxation 
dSCmax = 0;
for i=1:N
    dSCmax = max(abs(dU(2*i)),dSCmax);
end
dSCobj = 0.2;
alpha = min(1,dSCobj/dSCmax);
dU = alpha*dU;
alphamin = alpha;


//incrementation Newton 
for i=1:N 
 p(i) = p(i) + dU(2*i-1);
 if (DIPH(i)==0) then 
    C(i) = C(i) + dU(2*i);
    if (C(i) < f_cbar(p(i)) ) then
     DIPH(i) = 0;
     S(i) = 0;
    else
     DIPH(i) = 1;
     S(i) = 0;
     C(i) = f_cbar(p(i)); 
    end
    
 else
    S(i) = S(i) + dU(2*i) ;
    if (S(i) < 0) then
       DIPH(i)= 0 ;
       S(i) = 0;
       C(i) = f_cbar(p(i)) ;       
    else
       DIPH(i)= 1 ;
       C(i) = f_cbar(p(i)) ;       
    end
 end

end



//troncature :
for i=1:N
  //des saturations 
  if (S(i)>1)
    S(i) = 1;
  end
end 





//residu
[rhoo,rhog,dprhoo,dcrhoo,dprhog,dpcbar,fluxint,fluxwell,upwind,R] = residu(p,S,C,DIPH,p0,S0,C0,dt)

//norme du residu 
normR = 0;
for i=1:N
  normR = normR + R(i,1)**2 + R(i,2)**2;
end
normR = sqrt(normR);

yn(itn) = normR/normR0;
//yn(itn) = normR;

alphan(itn) = alphamin;

//fin Newton ---------------------
end
//fin Newton ---------------------

//restart pas de temps  
if (itn>Newtmax-1) 
 C = C0;
 S = S0;
 p = p0;
 DIPH = DIPH_0;
 t = t-dt;
 dt = dt/2;
 restart = 1
 n = n -1;
else 
 restart = 0;  
end  
 
//gestion pas de temps et calcul debits huile et gaz du well producteur
 if (restart==0) 
// debut if pas de restart ---------  
 temps(n) = temps(n-1) + dt;
 pasdetemps(n) = dt;
 C0=C;
 S0 =S;
 p0=p;
 rhoo0 = rhoo;
 rhog0 = rhog;
 DIPH_0 = DIPH ;
 dt = 1.2*dt;
 if (dt>dtmax)
  dt = dtmax;
 end
 if (t+dt>tf)
  dt = tf-t;
 end



 if (fluxwell(1)>0)
  m1 = cellwell(1) ;
  debito(n) =  fluxwell(1)*rhoo(m1)*Mo(1-S(m1))*(1-C(m1));  
  debitg(n) = fluxwell(1)*( rhoo(m1)*Mo(1-S(m1))*C(m1)+ rhog(m1)*Mg(S(m1)) );
 end
  debitTg(n) = debitTg(n-1) + debitg(n)*dt;  
  debitTo(n) = debitTo(n-1) + debito(n)*dt; 

//affichage saturation d'huile et pression 
if (int(n/1)*1==n)
   
scf(2)
//xbasc(2)
title('saturation en gaz') ;
plot(Xi,S)


scf(3)
//xbasc(3)
title('concentration gaz dans la phase huile') ;
plot(Xi,C)

scf(4)
//xbasc(4)
title('pression') ;
plot(Xi,p)


scf(9)
//xbasc(9)
title('indice de phase') ;
plot(Xi,DIPH)


//scf(10)
////xbasc(10)
//title('rhoo') ;
//plot(Xi,rhoo)
//
//
//scf(11)
////xbasc(11)
//title('rhog') ;
//plot(Xi,rhog)


end
//fin du if pour dessin

// 
end
// fin if pas de restart --------


//affichage courbes des residus Newton 
scf(26);
plot(xn,log10(yn),'b-*')
clear yn;

//scf(27);
//title('alphak')
//plot(xn,alphan,'b-*')
clear xn;
//clear alphan;

//fin boucle en temps ----------------------------------------
end
//fin boucle en temps ----------------------------------------

timer()



//affichages debits huile et gaz au producteur
scf(6);
title('debits leger (rouge) et lourd (vert) cumules') ;
plot(temps/3600/24,debitTg,'r-')
plot(temps/3600/24,debitTo,'g-')

scf(7);
title('debits leger (rouge) et lourd (vert) instantanné') ;
plot(temps/3600/24,debitg,'r-')
plot(temps/3600/24,debito,'g-')

//affichage marche en temps 
scf(8);
title('marche en temps') ;
plot(temps/3600/24,pasdetemps/3600/24,'-*')













