Why am I getting Nan values?

close all

clear all

dx=20;

imax=20;

dy=30;

dz=30;

dt=0.01;

tmax=10;

kfeff=100;

km=10;

phif=0.005;

phim=0.20;

GW=0.45;

GO=0.35;

l=3;

muo=0.8;

muw=0.5;

pi=2000;

pe=1600;

sigma=4*3/3^2;

sigmaz=4/3^2;

x=0.006328;

vi=dx*dy*dz;

vii=dx*dy*dz;

v=dx*dy*dz;

Swrf=0.05;

Sorf=0.05;

Swrm=0.35;

Sorm=0.3;

alpha=-0.55;

Swmx=0.55;

eps=0.0001;

KRWF=0.9;

KROF=0.9;

KRWM=0.07;

KROM=0.7;

nwf=2;

nwm=2;

nof=3;

nom=3;

qi=600;

%water in the fracture equations

Tfw=x*dy*dz/dx*kfeff/muw;

S=x*vi*sigma*km/muw;

%Oil in the fracture equations

Tfo=x*dy*dz/dx*kfeff/muo;

SO=x*vi*sigma*km/muo;

%intializing the system

for i=1:imax

Swf(i,1)=Swrf;

Swm(i,1)=Swrm;

pof(i,1)=1800;

pom(i,1)=1800;

end

for t=1:100

for i=1:imax

%oil saturation

Sof(i,t)=1-Swf(i,t);

Som(i,t)=1-Swm(i,t);

if Swf(i,t) <Swrf

Swf(i,t)=Swrf;

end

if Swm(i,t) <Swrm

Swm(i,t)=Swrm;

end

if Som(i,t) <Sorm

Som(i,t)=Sorm;

end

if Sof(i,t) <Sorf

Sof(i,t)=Sorf;

end

%hwf

hwf(i,t)=(Swf(i,t)-Swrf)/(1-Sorf-Swrf)*l;

%hwm

hwm(i,t)=(Swm(i,t)-Swrm)/(1-Sorm-Swrm)*l;

%Pcwom

pcwm(i,t)=2*log((1-Sorm-Swrm)/(Swm(i,t)-Swrm+eps));

%pcwof

pcwf(i,t)=0.1*log((1-Sorf-Swrf)/(Swf(i,t)-Swrf+eps));

%Relativie permeability in the fracture

%Krwf

krwf(i,t)=KRWF*((Swf(i,t)-Swrf)/(1-Sorf-Swrf))^nwf;

%Krof

krof(i,t)=KROF*((Sof(i,t)-Sorf)/(1-Sorf-Swrf))^nof;

%Relativie permeability in the matrix

%Krwm

krwm(i,t)=KRWM*((Swm(i,t)-Swrm)/(1-Sorm-Swrm))^nwm;

%Krom

krom(i,t)=KROM*((Som(i,t)-Sorm)/(1-Sorm-Swrm))^nom;

%calculating PI

WI=x*kfeff*krof(i,t)*muo*dy*dz*2/dx;

%relative permeability if statements

end

%%%%%%water in the fracture%%%%%

%FW

for i=1:imax-1

A(i,i+1)=Tfw*krwf(i,t);

end

%DW

for i=2:imax

A(i,i-1)=Tfw*krwf(i,t);

end

%EW

for i=1:imax

A(i,i)=-2*Tfw*krwf(i,t)-S*krwm(i,t);

end

A(1,1)=-Tfw*krwf(i,t)-S*krwm(i,t);

A(imax,imax)=-Tfw*krwf(i,t)-S*krwm(i,t);

%MFWR

for i=1:imax

A(i,i+20)= S*krwm(i,t);

end

%SWFR

for i=1:imax

A(i,i+40)= -vi*phif/dt;

end

%Rwf

for i=1:imax

if i==1

Rwf(i,1)=-Tfw*krwf(i,t)*pof(i+1,t)+Tfw*krwf(i,t)*pof(i,t)+Tfw*krwf(i,t)*pcwf(i+1,t)...

-Tfw*krwf(i,t)*pcwf(i,t)+Tfw*krwf(i,t)*pof(i,t)...

-Tfw*krwf(i,t)*pcwf(i,t)+S*krwm(i,t)*pof(i,t)...

-S*krwm(i,t)*pom(i,t)+S*krwm(i,t)*sigmaz*GW*hwf(i,t)/sigma-qi...

-S*krwm(i,t)*sigmaz*GW*hwm(i,t)/sigma-S*krwm(i,t)*pcwf(i,t)+S*krwm(i,t)*pcwm(i,t) ;

elseif i==imax

Rwf(i,1)=Tfw*krwf(i,t)*pof(i,t)...

-Tfw*krwf(i,t)*pcwf(i,t)+Tfw*krwf(i,t)*pof(i,t)-Tfw*krwf(i,t)*pof(i-1,t)...

-Tfw*krwf(i,t)*pcwf(i,t)+Tfw*krwf(i,t)*pcwf(i-1,t)+S*krwm(i,t)*pof(i,t)...

-S*krwm(i,t)*pom(i,t)+S*krwm(i,t)*sigmaz*GW*hwf(i,t)/sigma+WI*pe+WI*pof(i,t)...

-S*krwm(i,t)*sigmaz*GW*hwm(i,t)/sigma-S*krwm(i,t)*pcwf(i,t)+S*krwm(i,t)*pcwm(i,t);

else

Rwf(i,1)=-Tfw*krwf(i,t)*pof(i+1,t)+Tfw*krwf(i,t)*pof(i,t)+Tfw*krwf(i,t)*pcwf(i+1,t)...

-Tfw*krwf(i,t)*pcwf(i,t)+Tfw*krwf(i,t)*pof(i,t)-Tfw*krwf(i,t)*pof(i-1,t)...

-Tfw*krwf(i,t)*pcwf(i,t)+Tfw*krwf(i,t)*pcwf(i-1,t)+S*krwm(i,t)*pof(i,t)...

-S*krwm(i,t)*pom(i,t)+S*krwm(i,t)*sigmaz*GW*hwf(i,t)/sigma...

-S*krwm(i,t)*sigmaz*GW*hwm(i,t)/sigma-S*krwm(i,t)*pcwf(i,t)+S*krwm(i,t)*pcwm(i,t);

end

end

%%%%water in the matrix

for i=1:imax

%MW

A(i+20,i+20)=-x*sigma*km*krwm(i,t)*muw;

%SWMR

A(i+20,i+60)=-phim/dt;

%MFWL

A(i+20,i)=x*sigma*km*krwm(i,t)*muw;

end

%Rm

for i=1:imax

Rwm(i,1)=-x*sigma*km*krwm(i,t)*muw*sigmaz/sigma*GW*hwf(i,t)+x*sigma*km*krwm(i,t)*muw*sigmaz/sigma*GW*hwm(i,t)...

+x*pcwf(i,t)*sigma*km*krwm(i,t)*muw-x*sigma*km*krwm(i,t)*muw*pcwm(i,t)+x*sigma*km*krwm(i,t)*muw*pom(i,t)...

-x*sigma*km*krwm(i,t)*muw*pof(i,t);

end

%%%%%%Oil in the fracture%%%%%

%FO

for i=1:imax-1

A(i+40,i+1)=Tfo*krof(i,t);

end

%DO

for i=2:imax

A(i+40,i-1)=Tfo*krof(i,t);

end

%EO

for i=1:imax

A(i+40,i)=-2*Tfo*krof(i,t)-SO*krom(i,t);

end

A(1,1)=-Tfo*krof(i,t)-SO*krom(i,t);

A(imax,imax)=-Tfo*krof(i,t)-SO*krom(i,t);

%MFO

for i=1:imax

A(i+40,i+20)= SO*krom(i,t);

end

%SOF

for i=1:imax

A(i+40,i+40)= -vi*phif/dt;

end

%Rof

for i=1:imax

if i==1

Rof(1,1)=-Tfo*krof(i,t)*pof(i+1,t)+Tfo*krof(i,t)*pof(i,t)...

+Tfo*krof(i,t)*pof(i,t)...

+SO*krom(i,t)*pof(i,t)-SO*krom(i,t)*sigmaz*GO*hwm(i,t)/sigma...

-SO*krom(i,t)*pom(i,t)+SO*krom(i,t)*sigmaz*GO*hwf(i,t)/sigma...

-vi*phif/dt*Sof(i,t)+vi*phif/dt-vi*phif/dt*Swf(i,t)-qi;

elseif i==imax

Rof(imax,1)=Tfo*krof(i,t)*pof(i,t)...

+Tfo*krof(i,t)*pof(i,t)-Tfo*krof(i,t)*pof(i-1,t)...

+SO*krom(i,t)*pof(i,t)-SO*krom(i,t)*sigmaz*GO*hwm(i,t)/sigma...

-SO*krom(i,t)*pom(i,t)+SO*krom(i,t)*sigmaz*GO*hwf(i,t)/sigma...

-vi*phif/dt*Sof(i,t)+vi*phif/dt-vi*phif/dt*Swf(i,t)+WI*pe+WI*pof(i,t);

else

Rof(i,1)=-Tfo*krof(i,t)*pof(i+1,t)+Tfo*krof(i,t)*pof(i,t)...

+Tfo*krof(i,t)*pof(i,t)-Tfo*krof(i,t)*pof(i-1,t)...

+SO*krom(i,t)*pof(i,t)-SO*krom(i,t)*sigmaz*GO*hwm(i,t)/sigma...

-SO*krom(i,t)*pom(i,t)+SO*krom(i,t)*sigmaz*GO*hwf(i,t)/sigma...

-vi*phif/dt*Sof(i,t)+vi*phif/dt-vi*phif/dt*Swf(i,t);

end

end

%%%%%%Oil in the matrix%%%%%

%MO

for i=1:imax

A(i+60,i+20)=-x*sigma*km*krom(i,t)*muo;

%SOMR

A(i+60,i+60)=-phim/dt;

%MOFL

A(i+60,i)=x*sigma*km*krom(i,t)*muo;

end

%Rom

for i=1:imax

Rom(i,1)=-phim/dt*Som(i,t)-sigma*km*krom(i,t)*muo*sigmaz/sigma*GO*hwf(i,t)...

+sigma*km*krom(i,t)*muo*sigmaz/sigma*GO*hwm(i,t)...

-sigma*km*krom(i,t)*muo*pof(i,t)+sigma*km*krom(i,t)*muo*pom(i,t)...

+phim/dt-phim/dt*Swm(i,t);

end

%R matrix

for i=1:imax

R(i,1)=Rwf(i,1);

R(i+20,1)=Rwm(i,1);

R(i+40,1)=Rof(i,1);

R(i+60,1)=Rom(i,1);

end

delta=A\R;

for i=1:imax

pof(i,t+1)=pof(i,t)+delta(i,1);

pom(i,t+1)=pom(i,t)+delta(i+20,1);

Swf(i,t+1)=Swf(i,t)+delta(i+40,1);

Swm(i,t+1)=Swm(i,t)+delta(i+60,1);

Sof(i,t+1)=1-Swf(i,t+1);

Som(i,t+1)=1-Swm(i,t+1);

end

if Swf(i,t+1)<Swrf

krwf(i,t+1)=0;

elseif Swf(i,t+1)>1-Sorf

krwf(i,t+1)=KRWF;

end

if 1-Swf(i,t+1)<Sorf

krof(i,t+1)=0;

elseif 1-Swf(i,t+1)>1-Sorf

krof(i,t+1)=KROF;

end

if Swm(i,t+1)<Swrm

krwm(i,t+1)=0;

elseif Swf(i,t+1)>1-Sorf

krwm(i,t+1)=KRWM;

end

if 1-Swm(i,t+1)<Sorm

krom(i,t+1)=0;

elseif 1-Swm(i,t+1)>1-Sorm

krom(i,t+1)=KROM;

end

%hwf

hwf(i,t+1)=(Swf(i,t+1)-Swrf)/(1-Sorf-Swrf)*l;

%hwm

hwm(i,t+1)=(Swm(i,t+1)-Swrm)/(1-Sorm-Swrm)*l;

%Pcwom

pcwm(i,t+1)=2*log((1-Sorm-Swrm)/(Swm(i,t+1)-Swrm+eps));

%pcwof

pcwf(i,t+1)=0.1*log((1-Sorf-Swrf)/(Swf(i,t+1)-Swrf+eps));

%Relativie permeability in the fracture

%Krwf

krwf(i,t+1)=KRWF*((Swf(i,t+1)-Swrf)/(1-Sorf-Swrf))^nwf;

%Krof

krof(i,t+1)=KROF*(((1-Swf(i,t+1))-Sorf)/(1-Sorf-Swrf))^nof;

%Relativie permeability in the matrix

%Krwm

krwm(i,t+1)=KRWM*((Swm(i,t+1)-Swrm)/(1-Sorm-Swrm))^nwm;

%Krom

krom(i,t+1)=KROM*((1-(Swm(i,t+1))-Sorm)/(1-Sorm-Swrm))^nom;

%calculating PI

WI=x*kfeff*krof(i,t+1)*muo*dy*dz*2/dx;

end

[darova]

Can you be more specific? You have a long code

Please use attachment button for long code

Answer 1

[Constantino Carlos Reyes-Aldasoro]

Hello

There are several issues, but primarily about the way you use your code.

1) The code as it is is very long and convoluted, you could use functions and call these from your code and perhaps subfunctions.

2) Many of your for-loops could be replaced by proper matrix handling, that is if you want to multiply/add/subtract elements of one matrix by those of another matrix you do not need to loop over rows and columns, you can directly do A-B

Now, those two would make your code more efficient and clear, but not necessarily solve the problem at hand that is the NaNs.

After running your code and inspecting some of your variables (Som, Swm, ...) I can see that some of the variables have valid values (in fact the same) in rows 1:20 and columns 1:9 and then NaNs in columns 10:101. Thus the error must come in one of those loops, perhaps you are only running the first 9 columns and not the rest?

Now other variables like delta are just nans, same as Rom, Rwf, ...

My advice would be to reduce your problem, run for a small set, and calculate your variables one by one then you will be able to spot where the NaNs start

Answer 2

[Steven Lord]

Set an error breakpoint to stop when an Inf or NaN value gets generated then run your code. When MATLAB stops, look at the line that triggered the breakpoint. Most likely if that line generated an Inf your code overflowed or divided by zero. If it generated a NaN you probably computed 0/0.

Determine why the involved variables have the values that led to the nonfinite by setting breakpoints on earlier lines and rerunning your code, tracing back to where the code stops doing what you expected it to do.