Untitled

%%  Numerical Inversion using Gaver-Stehfest Algorithm [MATLAB]
%
% Advisor           : Prof. Asep K. Permadi, Ph.D
% Code Written by   : Billal Aslam, M.T.
% Modified By       : Kartika Fajarwati Hartono

%****************************************************************************
% Purpose           : Numerically inverse flow solution from Laplace Domain
%                     using Gvr-Sthfst Algorithm.
% Program Features  : Provide comparison with exact solution from 
%                     [0] Simple Exponential Functions  
%                     [1] Transient Bilinear flow in fractured well
%                     [2] Transient radial flow in no-flow bounded reservoir

%****************************************************************************

clc;
clearvars;
close all;

%% Main Calculation (Analytically Derived Solution/Numerical LT Inverse)

t = 1:1000:1000000; %time interval
n = 4; %stehfest sum (must be EVEN!) stevhest number
f=zeros(size(t,2),1);
G=zeros(size(t,2),1);
abse=zeros(size(t,2),1);


for i=1:size(t,2)
f(i) = feval('fun2',t(i));     %Analytical function 
G(i) = gavsteh('Lfun2',t(i),n); %Inverse Laplace function
abse(i)=abs(G(i)-f(i));        %Absolute error
end


%% Plot Setting
figure
% Cartesian plot (function #1 & #2)
 plot(t,f,'o'); hold on
 plot(t,G,'LineWidth',2);
 xlabel('t'); ylabel('f(t)');
 legend('analytic','GS approximation');
 figure
 plot(t,abse);

%log log plot (function #3)
loglog(t,G,'LineWidth',2);
xlabel('tD'); ylabel('pD(tD)');
title('pD(tD) for rD = 300')
%% Function Libraries
function f=fun1(t)
%simple function
f=exp(-t);
end

function f=Lfun1(s)
f=1/(s+1);
end

function f=fun2(t)
%Exact Solution of Pwf in bilinear transient flow fractured well
Pi  = 2200;
mu  = 12;
q   = 3200;
L   = 0.2;
H   = 2;
D   = 5e-8;
k   = 10;
C   = mu*q/(2*k*L*H);

f=Pi-C*sqrt(D*t/pi()); %Exact Solution equation
end

function f=Lfun2(s)
%Laplace Domain Solution of Pwf in bilinear transient flow fractured well
Pi  = 2200;
mu  = 12;
q   = 3200;
L   = 0.2;
H   = 2;
D   = 5e-8;
k   = 10;
C   = mu*q/(2*k*L*H);

f=(Pi/s)-(0.5*C*sqrt(D)*s^(-3/2)); %laplace equation 
end

function f=Lfun3(s)
%Laplace Domain Solution of PD in const. rate radial bounded reservoir
rD=1;
reD=1000;

f=(1/s^1.5)*((besselk(1,reD*sqrt(s))*besseli(0,rD*sqrt(s)))+...
    (besseli(1,reD*sqrt(s))*besselk(0,rD*sqrt(s))))/...
    (besseli(1,reD*sqrt(s))*besselk(1,sqrt(s))-...
    besselk(1,reD*sqrt(s))*besseli(1,sqrt(s))); %  bessel function
end


%% Laplace Transform Inversion Subroutine
%  ilt=gavsteh(funname,t,L)
%
%    funname       The name of the function to be transformed.
%    t             The transform argument (usually a snapshot of time).
%    ilt             The value of the inverse transform 
%    L              number of coefficient ---> depends on computer word length used
%                   (examples: L=8, 10, 12, 14, 16, so on..)
%
%  Wahyu Srigutomo
%  Physics Department, Bandung Institute of Tech., Indonesia, 2006
%  Numerical Inverse Laplace Transform using Gaver-Stehfest method
%
%References:
% 1. Stehfest, H., 1970, Algorithm 368: Numerical inversion of Laplace transform,
%    Communication of the ACM, vol. 13 no. 1 p. 47-49

function ilt=gavsteh(funname,t,L)
nn2 = L/2;
%nn21= nn2+1;

for n = 1:L
    z = 0.0;
	for k = floor( ( n + 1 ) / 2 ):min(n,nn2)
         z = z + ((k^nn2)*factorial(2*k))/ ...
            (factorial(nn2-k)*factorial(k)*factorial(k-1)* ...
            factorial(n-k)*factorial(2*k - n));
    end 
    v(n)=(-1)^(n+nn2)*z;
end

sum = 0.0;
ln2_on_t = log(2.0) / t;
for n = 1:L
    p = n * ln2_on_t;
    sum = sum + v(n) * feval(funname,p);
end 
	ilt = sum * ln2_on_t;
end