
%%% this function solves the couled partial differential equations of
%%% subsurface tracer transport
function [sim_t,Ca,Cb,Cc]=SolvePDE(w, Q, q, constant)

Q_inj1=Q(1);
Q_inj2=Q(2);
Q_pro=-Q(3);
q1=q(1)*w(4);
q2=q(2)*w(5);

T1=q1/Q_inj1;
T2=q2/Q_inj2;


K=constant(1);
h=constant(2);
phi=constant(3);
T_s=constant(4);
CA_0=constant(5);


Sor=w(1);
k=w(2);
alpha=w(3);

    
beta=K*Sor/(1-Sor);




m = 1;  %%%% define geometry number see pdepe referece page for details
x = linspace(0.5,100,200); 
xout=x;

phase=1;

t = linspace(0,T1,3);  %%%% define time points to calculate concentration
sol = pdepe(m,@allpde,@ic,@bc,x,t);   %%%% define the pde, initial,
output1=sol(end,:,1);
output2=sol(end,:,2);
output3=sol(end,:,3);

phase=2;
t = linspace(0,T2,3);  %%%% define time points to calculate concentration
sol2 = pdepe(m,@allpde,@ic,@bc,x,t);   %%%% define the pde, initial,
output1=sol2(end,:,1);
output2=sol2(end,:,2);
output3=sol2(end,:,3);

phase=3;
t = linspace(0,T_s,3);  %%%% define time points to calculate concentration
sol3 = pdepe(m,@allpde,@ic,@bc,x,t);   %%%% define the pde, initial,
output1=sol3(end,:,1);
output2=sol3(end,:,2);
output3=sol3(end,:,3);
% 
phase=4;
t = linspace(0,2.5,floor(2.5*50));  %%%% define time points to calculate concentration
sol4 = pdepe(m,@allpde,@ic,@bc,x,t);   %%%% define the pde, initial,


sim_t=-t*Q_pro;



Ca=CA_0*sol4(:,1,1);
Cb=CA_0*sol4(:,1,2);
Cc=CA_0*sol4(:,1,3);

function [c,f,s] = allpde(x,t,u,DuDx)
    
 
switch(phase)   
    %%%%% different operation stages, two steps of injection, well shut-in,
    %%%%% and production
    case 1
        Q0=Q_inj1;
    case 2
        Q0=Q_inj2;
    case 3
        Q0=0;      
    case 4
        Q0=Q_pro;
end
vr=5.615*Q0/(2*pi*x*h*phi*(1-Sor));
Kr=alpha*abs(vr);
c = [1+beta;1;1];
f = [Kr*DuDx(1);Kr*DuDx(2);Kr*DuDx(3)];
s = [-k*u(1)-vr*DuDx(1);k*u(1)-vr*DuDx(2);-vr*DuDx(3)];
end

function u0 = ic(x)
switch(phase)
    case 1
        u0 = [0;0;0];
    otherwise
        u0=[output1(xout==x);output2(xout==x);output3(xout==x)] ;
end
end


function [pl,ql,pr,qr] = bc(xl,ul,xr,ur,t)
switch(phase)
    case 1
        pl=[ul(1)-1;ul(2)-0;ul(3)-1];
        ql=[0;0;0];
    case 2
        pl=[ul(1);ul(2);ul(3)-1];
        ql=[0;0;0];
    otherwise
        pl=[0;0;0];
        ql=[1;1;1];
end
 pr=[0;0;0];
 qr=[1;1;1];

end






end




  