get_start.m

Petroleum Reservoir Simulation

%Fill in every line and get pwi, swi, phi and so on

%function [JM, fwi,fni]=get_start(pwi,swi,phi,pwit,swit,phit,nr)

for ii=1:nr
   pci(ii)=aa.*exp(-(swi(ii)./bb));
   pci1(ii)=-pci(ii)./bb;
   pci11(ii)=pci(ii)./bb./bb;
end

for ii=1:nr
   if ii<nr
      pcip1(ii)=pci(ii+1);
   end
   pcip1(nr)=0;
   
   pcim1(1)=0;
   pcim11(1)=0;
   pcim111(1)=0;
   if ii>1
      pcim1(ii)=pci(ii-1);
      pcim11(ii)=-pcim1(ii)./bb;
      pcim111(ii)=pcim1(ii)./bb./bb;
   end
end

for ii=1:nr
   if ii<nr
      pwip1(ii)=pwi(ii+1);
      swip1(ii)=swi(ii+1);
   end
   pwip1(nr)=0;
   swip1(nr)=0;
   
   pwim1(1)=0;
   swim1(1)=0;
   if ii>1
      pwim1(ii)=pwi(ii-1);
      swim1(ii)=swi(ii-1);
   end
end

pni=pci+pwi;

   

%fill in "TW"   
%twip[   ]   Twi+     
for ii=1:nr
   krwi(ii)=swi(ii).^4;  %get "swi" first
   krwi1(ii)=4.*(swi(ii).^3);
   krni(ii)=(1-swi(ii)).^2;  %get "swi" first
   krni1(ii)=-2.*(1-swi(ii));
   
   bwi(ii)=bwr+bw1.*(pwi(ii)-pwr);   %get "pwi" first 
   twip(ii)=kw(ii).*krwi(ii).*bwi(ii);
   bni(ii)=bnr+bn1.*(pni(ii)-pnr);   %get "pwi" first 
   tnip(ii)=kn(ii).*krni(ii).*bni(ii);
   
   %phi(ii)=phir.*exp(cf.*(pwi(ii)-pr));
   phitp(ii)=phir.*exp(cf.*(pwi(ii)-pr));
   phi1(ii)=cf.*phir.*exp(cf.*(pwi(ii)-pr));
   phi11(ii)=cf.*cf.*phir.*exp(cf.*(pwi(ii)-pr));

if ii==nr
   twip(ii)=0;
   tnip(ii)=0;
end


end

%twim[   ]   Twi-   krwim1[   ]   Krwi-1     bwim1[   ]   bwi-1
for ii=1:nr
   
 if ii==1
    twim(ii)=0;
    krwim1(ii)=0;
    krwim11(ii)=0;
    bwim1(ii)=0;
    
    tnim(ii)=0;
    krnim1(ii)=0;
    krnim11(ii)=0;
    bnim1(ii)=0;
 else
    twim(ii)=kw(ii).*krwi(ii-1).*bwi(ii-1);
    krwim1(ii)=krwi(ii-1);
    krwim11(ii)=krwi1(ii-1);
    bwim1(ii)=bwi(ii-1);
    
    tnim(ii)=kn(ii).*krni(ii-1).*bni(ii-1);
    krnim1(ii)=krni(ii-1);
    krnim11(ii)=krni1(ii-1);
    bnim1(ii)=bni(ii-1);
end

end

bwitp=bwi;
bnitp=bni;
phitp=phi;



fwi1=twip.*pwip1-(twip+twim+vb./dt.*(swit.*phit).*bwi1+vb./dt.*(swit.*bwitp.*phi1)).*pwi;
fwi2=twim.*pwim1-vb./dt.*(bwitp.*phitp).*swi+vb./dt.*(bwitp.*phitp).*swit-twim.*rw.*(zi-zim1)+qwi;
fwi3=(vb./dt.*(swit.*phit).*bwi1+vb./dt.*(swit.*bwitp.*phi1)).*pwit+twip.*rw.*(zip1-zi);
fwi=fwi1+fwi2+fwi3;

fwp1=-kw.*krwim1.*bwim11.*pwi+kw.*krwim1.*bwim11.*pwim1+kw.*krwim1.*bwim1-kw.*krwim1.*bwim11.*rw.*(zi-zim1);
fws1=kw.*krwim11.*bwim1.*(pwim1-pwi-rw.*(zi-zim1));

fwp21=kw.*krwi.*bwi1.*pwip1-kw.*krwi.*bwi1.*pwi-kw.*krwi.*bwitp-(kw.*krwim1.*bwim1+vb./dt.*swit.*phit.*bwi1);
fwp22=-vb./dt.*swit.*(bwi1.*phi1.*pwi+bwitp.*phi11.*pwi+bwitp.*phi1)-vb./dt.*swi.*(bwitp.*phi1+bwi1.*phi);
fwp23=vb./dt.*swit.*(bwitp.*phi1+bwi1.*phi)+vb./dt.*pwit.*swit.*(bwi1.*phi1+bwitp.*phi11)+kw.*krwi.*bwi1.*rw.*(zip1-zi);
fwp2=fwp21+fwp22+fwp23;
fws2=kw.*krwi1.*bwi.*(pwip1-pwi+rw.*(zip1-zi))-vb./dt.*(bwitp.*phitp);

fwp3=kw.*krwi.*bwi;
for ii=1:nr
   fws3(ii)=0;
end




fni1=tnip.*pwip1-(tnip+tnim+vb./dt.*((1-swit).*phit).*bni1+vb./dt.*((1-swit).*bnitp.*phi1)).*pwi;
fni2=tnim.*pwim1+tnip.*pci1.*swip1+tnim.*pcim11.*swim1-tnim.*rn.*(zi-zim1)+qni;
fni3=-(tnip.*pci1+tnim.*pcim11+vb./dt.*((1-swit).*phit).*bni1.*pci1-vb./dt.*(phitp.*bnitp)).*swi;
fni4=(vb./dt.*((1-swit).*phit).*bni1+vb./dt.*((1-swit).*bnitp.*phi1)).*pwit;
fni5=(vb./dt.*((1-swit).*phit).*bni1.*pci1-vb./dt.*(phitp.*bnitp)).*swit+tnip.*rn.*(zip1-zi);
fni=fni1+fni2+fni3+fni4+fni5;

fnp11=-kn.*krnim1.*bnim11.*pwi+kn.*krnim1.*bnim11.*pwim1+kn.*krnim1.*bnim1;
fnp12=kn.*krnim1.*bnim11.*pcim11.*swim1-kn.*krnim1.*bnim11.*pcim11.*swi-kn.*krnim1.*bnim11.*rn.*(zi-zim1);
fnp1=fnp11+fnp12;
fns11=kn.*krnim11.*bnim1.*(pwim1-pwi-rn.*(zi-zim1))+kn.*krnim11.*bnim1.*pcim11.*swim1;
fns12=kn.*krnim1.*bnim1.*pcim111.*swim1+kn.*krnim1.*bnim1.*pcim11-kn.*krnim11.*bnim1.*pcim11.*swi-kn.*krnim1.*bnim1.*pcim111.*swi;
fns1=fns11+fns12;

fnp21=kn.*krni.*bni1.*pwip1-kn.*krni.*bni1.*pwi-kn.*krni.*bni-(kn.*krnim1.*bnim1+vb./dt.*(1-swit).*phit.*bni1);
fnp22=-vb./dt.*(1-swit).*(bni1.*phi1.*pwi+bnitp.*phi11.*pwi+bnitp.*phi1)+kn.*krni.*bni1.*pci1.*swip1;
fnp23=-kn.*krni.*bni1.*pci1.*swi+vb./dt.*swi.*(bnitp.*phi1+bni1.*phi)+vb./dt.*(1-swit).*pwit.*(bni1.*phi1+bnitp.*phi11);
fnp24=-vb./dt.*swit.*(bnitp.*phi1+bni1.*phi)+kn.*krni.*bni1.*rn.*(zip1-zi);
fnp2=fnp21+fnp22+fnp23+fnp24;
fns21=kn.*krni1.*bni.*(pwip1-pwi+rn.*(zip1-zi))+swip1.*(kn.*krni1.*bni.*pci1+kn.*krni.*bni.*pci11);
fns22=-kn.*krni1.*bni.*pci1.*swi-kn.*krni.*bni.*swi.*pci11-kn.*krni.*bni.*pci1+vb./dt.*((1-swit).*phit).*bni1.*swit.*pci11;
fns23=-vb./dt.*((1-swit).*phit).*bni1.*(pci11.*swi+pci1)-(kn.*krnim1.*bnim1.*pcim11-vb./dt.*(phitp.*bnitp));
fns2=fns21+fns22+fns23;

fnp3=kn.*krni.*bni;
fns3=kn.*krni.*bni.*pci1;

fwp=flx_fill_in(fwp1,fwp2,fwp3,nr);
fws=flx_fill_in(fws1,fws2,fws3,nr);
fw=fl_fill_in(fwp,fws,nr);

fnp=flx_fill_in(fnp1,fnp2,fnp3,nr);
fns=flx_fill_in(fns1,fns2,fns3,nr);
fn=fl_fill_in(fnp,fns,nr);

JM=jacobi(fw,fn,nr); %get the jacobi matrix
%JM(2*nr-3,2*nr-1)=1000;
%JM(2*nr-2,2*nr-1)=1000;
%JM(2*nr-1,2*nr-1)=1000;
%JM(2*nr,2*nr-1)=1000;



fk=fl_fill_in(fwi,fni,nr); %get the right Fk
%fk(2*nr-1)=1000.*(pwi(nr)-pwf);
%fk(2*nr)=1000.*(pni(nr)-pwf);

fki=fk;