close all
clear all

% This matlab code calculates the upscaled flux F and integrated saturation
% S for the case of sedimentary strata, then uses these to solve the
% similarity equation to give the gravity current thickness h, as well as
% the prefactor etaN

Bo=1;       % Bond number
Hrat=1;     % Width ratio of layers
krat=0.1;   % Permeability ratio
lambda=3;   % Capillary pressure exponent
alpha=4;    % Relative permeability exponent 
b=0.43;     % Pore entry pressure-permeability power law
a=0.14;     % Porosity-permeability power law

fprintf('Building flux...\n')

Krnfun=@(s) s.^alpha ;

n=1e3;                  % # Points in 'horizontal'
m=1e3;                  % # Points in vertical
hmax=5*(1/sqrt(Bo)+1);  % Approximate scaling for h range

h=linspace(0,hmax,n);

if Hrat<1
    Hrat=1/Hrat;
    krat=1/krat;
end

% Calculate pore entry pressure, permeability and porosity within layers

A1=Hrat/(1+Hrat);
A2=krat;

Peval1=(1/(A1+(1-A1)*A2))^(-b);
Peval2=(A2/(A1+(1-A1)*A2))^(-b);

K1=(1/(A1+(1-A1)*A2));
K2=(A2/(A1+(1-A1)*A2));

Phi1=(1/(A1+(1-A1)*A2))^a;
Phi2=(A2/(A1+(1-A1)*A2))^a;

Phitot=A1*Phi1+(1-A1)*Phi2;
Phi1=Phi1/Phitot;
Phi2=Phi2/Phitot;

Ktot=A1*K1+(1-A1)*K2;
K1=K1/Ktot;
K2=K2/Ktot;

for i=1:n
    
    if A2<1
        H=h(i);
        Z=flip(H-logspace(-6,log10(H),m));
        Z(1)=0;
        
        H2=max(H-1/Bo*(Peval2-Peval1),0);
        
        Z2=flip(H2-logspace(-6,log10(H2),m));
        Z2(1)=0;
        
        Sc1=max(1-(Peval1./Peval1 + Bo*(H-Z)./Peval1).^(-lambda),0);
        DScDh1=lambda*Bo./Peval1.*(Peval1./Peval1 + Bo*(H-Z)./Peval1).^(-lambda-1);
        DScDh1(Sc1==0)=0;
        
        Sc2=max(1-(Peval1./Peval2 + Bo*(H-Z2)./Peval2).^(-lambda),0);
        DScDh2=lambda*Bo./Peval2.*(Peval1./Peval2 + Bo*(H-Z2)./Peval2).^(-lambda-1);
        DScDh2(Sc2==0)=0;
        
        Scbar1(i)=1/H*trapz(Z,Sc1.*Phi1);
        Scdashbar1(i)=1/H*trapz(Z,DScDh1.*Phi1);
        krnc1(i)=1/H*trapz(Z,K1.*Krnfun(Sc1));
        
        Scbar2(i)=1/H*trapz(Z2,Sc2.*Phi2);
        Scdashbar2(i)=1/H*trapz(Z2,DScDh2.*Phi2);
        krnc2(i)=1/H*trapz(Z2,K2.*Krnfun(Sc2));
    else
        H=h(i);
        Z=flip(H-logspace(-6,log10(H),m));
        Z(1)=0;
        
        H2=max(H-1/Bo*(Peval1-Peval2),0);
        
        Z2=flip(H2-logspace(-6,log10(H2),m));
        Z2(1)=0;
        
        Sc1=max(1-(Peval2./Peval1 + Bo*(H-Z2)./Peval1).^(-lambda),0);
        DScDh1=lambda*Bo./Peval1.*(Peval2./Peval1 + Bo*(H-Z2)./Peval1).^(-lambda-1);
        DScDh1(Sc1==0)=0;
        
        Sc2=max(1-(Peval2./Peval2 + Bo*(H-Z)./Peval2).^(-lambda),0);
        DScDh2=lambda*Bo./Peval2.*(Peval2./Peval2 + Bo*(H-Z)./Peval2).^(-lambda-1);
        DScDh2(Sc2==0)=0;
        
        Scbar1(i)=1/H*trapz(Z2,Sc1.*Phi1);
        Scdashbar1(i)=1/H*trapz(Z2,DScDh1.*Phi1);
        krnc1(i)=1/H*trapz(Z2,K1.*Krnfun(Sc1));
        
        Scbar2(i)=1/H*trapz(Z,Sc2.*Phi2);
        Scdashbar2(i)=1/H*trapz(Z,DScDh2.*Phi2);
        krnc2(i)=1/H*trapz(Z,K2.*Krnfun(Sc2));
        
    end
    
end

K=A1*krnc1+(1-A1)*krnc2;
meansat=A1*Scbar1+(1-A1)*Scbar2;
meansatdash=A1*Scdashbar1+(1-A1)*Scdashbar2;

Flux=K./meansatdash;
Sat=meansat.*h;

Flux(1)=0;
Sat(1)=0;

save('flux','Sat','Flux','h','meansatdash');

fprintf('Caclulating prefactor...\n')

f=@(etaN) etafind(etaN,Bo);
etaN=fsolve(f,1);

[out,hsol,Sinv] = etafind(etaN,Bo);
N=length(hsol);

hmid=hsol(N/2);
hsol=hsol/(2*hmid);

save('solution','hmid','etaN','hsol')

figure(1);clf;
y=linspace(0,1,N);
plot(y,hsol,'k','linewidth',2)
grid on
box on
set(gca,'fontsize',20)
xlabel('$$y$$','interpreter','latex')
ylabel('$$h$$','interpreter','latex')
yax=ylim;
ylim([0,yax(2)])
xlim([0,1])

figure(2);clf;
y=linspace(0,1,N);
plot(h(h<max(hsol)),Sat(h<max(hsol)),'k','linewidth',2)
grid on
box on
set(gca,'fontsize',20)
xlabel('$$h$$','interpreter','latex')
ylabel('$$\mathcal{S}$$','interpreter','latex')

figure(3);clf;
y=linspace(0,1,N);
plot(Sat(h<max(hsol)),Flux(h<max(hsol)),'k','linewidth',2)
grid on
box on
set(gca,'fontsize',20)
xlabel('$$\mathcal{S}$$','interpreter','latex')
ylabel('$$\mathcal{F}$$','interpreter','latex')


fprintf(strcat('\neta_N=',num2str(etaN),'\n'))
fprintf(strcat('\nh_half=',num2str(hmid),'\n'))

function [out,h,Sinv] = etafind(etaN,Bo)

load('flux');

Sinv=@(S) interp1(Sat,h,S,'linear','extrap');
flux=@(f) interp1(Sat,Flux,f,'linear','extrap');

N=10000;
L=-0.001;
f=(1/sqrt(Bo)+1)/1000;
yi=1;

rhsf=L/(flux(f)*yi);
rhsL=-1/2*etaN^2*yi^2*L/(flux(f)*yi);

y=linspace(0,1,N);
dy=y(2)-y(1);
fs=[f];
for i=1:N-1
    yi=y(N-i);
    
    f=f-dy*rhsf;
    L=L-dy*rhsL;
    
    rhsf=L/(flux(f)*yi);
    rhsL=-1/2*etaN^2*yi^2*L/(flux(f)*yi);
    
    fs=[fs,f];
end
fs=flip(fs);
intval=etaN^2*2*pi*trapz(y,y.*fs);
out=intval-1;
h=Sinv(fs);
end
