Error when using erfc. Input must be real and full...

50 visualizaciones (últimos 30 días)
Svante Monie
Svante Monie el 17 de Ag. de 2015
Comentada: Svante Monie el 20 de Ag. de 2015
I really would appreciate some help here...
In the attached code I have a problem with a variable, s, which is only valid as 0<s<pi. I tried to make it as a vector, but the the stepsize seemed to cause problems. When I worked the code against a reference case it seemed like if calculations above a time span of 100 years tolerated the same stepsize, but every year less than 100 years had more or less different tolerance in stepsize... Very disturbing. After tedious work I'd defined a vector for all years between 1 and 100 years.
So, when I then ran the code with actual data from field I got the same problem again: Error when using erfc. Input must be real and full. Even though I tried to have a longer timespan the problem persisted.
I've tried also to define the variable, s, like this:
s = [0.000001 pi]
But the problem persist.
I'm out of ideas now, anyone?
mvh
Svante W Monie
% !!Temporary parameter input for debugging!!
t_TOT = 100;
L = 933;
D = 750;
H = 200;
T = 320;
P = 161;
Lambda_c = 2.2;
Lambda_b = 2;
C_c = 2.2*10^6;
C_b = 2.295*10^6;
Q_w = 0.05141;
Por = 0.05;
C_s = 6.16*10^6;
Delta_T = 240;
C = 2.036175*10^6;
C_w = 4.1587*10^6;
% Thermal radius of one years reinjection, R_ty, and total radius of simulated time span, R_T.
t_y = 31536000; % t_y = one year in seconds
V_wy = Q_w*t_y; % [m^3] in one year.
R_Ty =sqrt((V_wy*C_w)/(pi * H * C)); % (8.1.31) p 8.16 Claesson et al
R_T = R_Ty*sqrt(t_TOT/1);
% Thermal shift in one year due to regional ground water
L_Ty = C_w*q_w/C; % (9.2.7) p. 9.4 Claesson et al
if L_Ty/R_Ty < 0.5; % (8.1.35) Claesson et al
L_Ty = 0;
q_wp = V_wy/(2 * pi * H * L * t_y); % Water flow in aquifer due to pump flow (when regional ground water flow can be neglected).
disp (' ');
disp ('The process is dominated by the reinjected water, i.e. the regional ground water');
disp (['flow can be neglected. Water flow will be determined of the pump flow: ', num2str(q_wp) ' [m^3/s]']);
else if L_Ty/R_Ty > 2; % (8.1.35) Claesson et al
disp ('Process is dominated by regional ground water flow.');
disp ('Simulation might not be valid!');
end
if 0.5 < L_Ty/R_Ty < 2; % (8.1.35) Claesson et al
disp ('The water flow process can be dominated by regional ground water flow. More');
disp ('detailed investigation is recommended.');
end
end
if R_Ty/L >= 0.75;
disp (' ');
disp ('Distance between wells is too short for specified flow rate.');
disp ('Thermal radius ~ distance between wells');
disp ('There is a risk of thermal short circuit!');
else
disp (' ');
disp ('Flow rate of injected water is not running risk of thermal short circuiting.');
disp ('Thermal radius << distance between wells');
end
disp (' ');
disp (['Thermal radius of one years injection will be ', num2str(R_Ty) ' [m].']);
disp (' ');
disp (['Thermal radius for entire time span, ', num2str(t_TOT) ' [year], will be ', num2str(R_T) ' [m]']);
% Defining time constraint for independence of influence from surface.
t_bp = (((2 * D) + H)^2/(pi *(Lambda_c/C_c)))/t_y;
% Defining melting time for thermal fluid in years.
t_m = (H^2*(2*C/(sqrt(Lambda_c*C_c)+sqrt(Lambda_b*C_b)))^2)/t_y; % (9.5.2.5) p. 9.42 Claesson et al
% Defining break-through time for injected cold water front at production well in years.
t_bt = ((pi*H*C*L^2)/(3*Q_w*C_w))/t_y; % [years] (9.6.1.2) p. 9.47 Claesson et al
disp (' ');
disp ('Break through time for first fraction of the cold water front to reach');
disp (['the production well will be ', num2str(t_bt) ' [year] and the time for']);
disp (['melting will be ', num2str(t_m) ' [year].']);
% Do computation of dimensionless factor gamma:
gamma = sqrt(t_bt/t_m); % (9.6.1.6) p. 9.48 Claesson et al
% Thermal influence radius
R_f = ((Q_w*C_w)/(2*pi*(Lambda_c+Lambda_b))); % (9.5.1.2) p. 9.34 Claesson et al
e = 0;
while e == 0;
% Definition of the thermal fluids circular paths
% if t_TOT > 100;
% S = 0.01;
% end
% if t_TOT <= 100;
% s_index = [0.3500 0.2650 0.1209 0.0350 0.0300 0.0310 0.0300 0.0220 0.0210 0.020 0.0210 0.0190 0.0180 0.0160 0.0140 0.0140 0.0140 0.0140 0.0140 0.0140 0.0130 0.0130 0.0130 0.0130 0.0130 0.0120 0.0120 0.0120 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0100 0.0100 0.0070 0.0040 0.0040 0.0040 0.0040 0.0040 0.0040 0.0030 0.0020 0.0100 0.0100 0.0100 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0120 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0110 0.0100 0.0100 0.0100 0.0100 0.0100 0.0100];
% S = s_index(t_TOT);
% end
%
% s = 0:S:pi;
S=0:0.01:pi;
s=zeros();
for i=2:length(S);
s(i-1)=S(i);
end;
%s = 0:0.1:pi;
% Definition of length for Tao vector
j = t_TOT/(length(s)-1);
t = (0:j:t_TOT);
% Do computation for dimensionless time Tao:
Tao = zeros();
for i=1:length(t);
Tao(i)=t(i)/t_bt;
end
% Determine when Tao >= 1
n=0;
u = 0;
for i=1:length(Tao);
n = Tao(i);
u = u + 1;
if n >= 1 || u==length(Tao);
break
end
end
if u == length(Tao);
disp(' ');
disp ('Time span is too short for disturbance to be detected. Try again with longer');
disp (['time span, not longer than the time restain of ' num2str(t_bp) 'years']);
disp ('after which impact of surface must be considered.');
% Create a dialog box for the user to input data.
prompt = {['Enter new time span for model, larger than ' num2str(t_TOT) ' years, less than ' num2str(t_bp) ':']};
Title = 'New time span for modeling';
answer = inputdlg(prompt,Title);
t_TOT = str2double(answer{1});
e = 0;
end
if u < length(Tao);
e = 1;
end
end
% Define f0(x) as a function for the throughput time between wells:
f0 = @(s) 3.*((sin(s)-(s.*cos(s)))./(sin(s)).^3);
% Calculation of temperature disturbance at production well
f_inv = zeros();
U_out = zeros();
for i=u:length(Tao);
c = Tao(i);
% Define fTao(s,c) as a function for the dimensionless time integral limit
% Note: s is only defined between 0 and pi.
FTao = @(s,c) (3.*((sin(s)-(s.*cos(s)))./(sin(s)).^3))-c;
fTao = @(s) FTao(s,c);
x0 = [0.00000000001 3.1415];
f_inv(i) = fzero(fTao,x0);
% Define the function f using f0(x) from above:
f = @(s) erfc(gamma*(f0(s)./sqrt(c-f0(s))));
% Compute the integral (note, the right boundary should be the inverse function of Tao):
U_out(i) = 1/pi * integral(f, 0, f_inv(i));
end
Temp_decline = Delta_T*U_out;
Year = Tao.*t_bt;
figure(1)
plot(Tao,U_out);
title ('Temperature disturbance at production well',...
'Fontsize',14);
xlabel ('Dimensionless time [\it\tau = t / t_{bt}\rm]','Fontsize',12);
ylabel('U_{out}(\tau)','Fontsize',12);
set(gca, 'Fontsize',10,...
'TickDir','out',...
'TickLength',[.02,.02],...
'YLim',[0,1])
grid on;
figure(2)
plot(Year,-Temp_decline);
hleg=legend(['With \DeltaT, between wells: ', num2str(Delta_T) ' °C'],...
'Location','NorthEast');
set(hleg,'FontAngle','italic','TextColor',[.3,.2,.1]);
title ('Temperature decline in reservoir',...
'Fontsize',14);
xlabel ('Years after comissioning','Fontsize',12);
ylabel('Degrees celsius [°C]','Fontsize',12);
set(gca, 'Fontsize',10,...
'TickDir','out',...
'TickLength',[.02,.02])
grid on;

Respuesta aceptada

Satyajeet Sasmal
Satyajeet Sasmal el 19 de Ag. de 2015
Hi Svante,
In your code, when you calculate the square root of a number, you need to make sure that he number is a real and positive number (as sqrt of negative numbers will give a complex number). Also, make sure that when you divide this sqrt value, make sure that "Division by Zero" does not happen. In such a case, set the value of c-f0(s) = 1 .
This solved the issue. I have attached the code snippet of where you need to bring in the required changes.
f_inv(i) = fzero(fTao,x0);
% Define the function f using f0(x) from above:
if (c-f0(s) > 0 ) % Check to see if this number is real and positive before calculating sqrt
con = c-f0(s);
else
con = 1;
end
f = @(s) erfc(gamma*(f0(s)./sqrt(con)));
Hope this was helpful!
  1 comentario
Svante Monie
Svante Monie el 20 de Ag. de 2015
Hi Satajeet
You are completely correct in your observation regarding that I must ensure that the input is positive and also that I can't get negative numbers under the square root in order to avoid complex numbers. In order to avoid this there is a section in the code (a massive an not so neat script, I know...) where this has been taken care of.
When I checked the code step-by-step I could confirm that this did work, so no negative and complex numbers where fed in to erfc...
Instead I discovered that if I did break it down as follows it worked just fine (discovered this just last night...) and so there was no issue with the parameter 's' as I first thought!
Many thanks for you comment, though - it's always nice to get response from all over the globe :-)
% Define the function f using f0(x) from above:
z = gamma(k)*(f0(s(i-u(k)+1))./sqrt(c-f0(s(i-u(k)+1))));
f = @(z) erfc(z);

Iniciar sesión para comentar.

Más respuestas (1)

Steven Lord
Steven Lord el 19 de Ag. de 2015
I didn't read through the whole code, but I did notice this:
if 0.5 < L_Ty/R_Ty < 2;
That doesn't do what you think it does. That condition is always true regardless of the value of L_ty and R_ty. Break it down into pieces:
if ((0.5 < L_Ty/R_Ty) < 2);
The expression in the innermost parentheses is either true (1) or false (0). Both 1 and 0 are less than 2. In order for that condition to be satisfied only when L_Ty/R_Ty is between 0.5 and 2, you need:
if (0.5 < L_Ty/R_Ty) & (L_Ty/R_Ty < 2);
I didn't check if L_Ty and R_Ty were scalars. If they are not, you would probably want to wrap those conditions in ANY or ALL depending on whether you want the body of the IF to be executed if any of the elements are between 0.5 and 2 or all are between 0.5 and 2.

Categorías

Más información sobre Data Type Identification en Help Center y File Exchange.

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by