Why am I getting Nan values?

12 views (last 30 days)
Sarah A Alruwayi
Sarah A Alruwayi on 1 Apr 2020
Answered: Steven Lord on 1 Apr 2020
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
  1 Comment
darova
darova on 1 Apr 2020
Can you be more specific? You have a long code
Please use attachment button for long code

Sign in to comment.

Answers (2)

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

Steven Lord
Steven Lord on 1 Apr 2020
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!