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

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;

Answer 1

[Satyajeet Sasmal]

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!

[Svante Monie]

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);

Answer 2

[Steven Lord]

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.

[Svante Monie]

Thank you! A relevant comment! I adopt it right away :-)