
%
%   Code to compute random permeability field for single-phase flow problem.
%

  clear all
  close all

%
% Set length scales and compute covariance matrix and weights.
%

  Lx = 0.5;
  Ly = 0.25;

  Nmodes = 104;
  Nquadx = 50;
  Nquady = 20;
  N = Nquadx*Nquady;
  hx = 2/Nquadx;
  hy = 1/Nquady;
  wx = hx*ones(1,Nquadx);
  wy = hy*ones(1,Nquady);
  w = hx*hy*ones(1,Nquadx*Nquady);  

  for j=1:Nquadx
    valuesx(j) = 2*j-1;
    x(j) = -1 + (hx/2)*valuesx(j);
  end
  for j=1:Nquady
    valuesy(j) = 2*j-1;
    y(j) = (hy/2)*valuesy(j);
  end
   
  for i=1:Nquadx
    for j=1:Nquadx
      Cx(i,j) = exp(-abs(x(j) - x(i))/Lx);
    end
  end

  for i=1:Nquady
    for j=1:Nquady
      Cy(i,j) = exp(-abs(y(j) - y(i))/Ly);
    end
  end

  indx = 1;
  for i=1:Nquadx
    for j=1:Nquady
      indy = 1;
      for k = 1:Nquadx
        for m = 1:Nquady
          C(indy,indx) = exp(-(abs(x(i) - x(k))/Lx) - (abs(y(j) - y(m))/Ly));
          indy = indy+1;
        end
      end
      indx = indx+1;
    end
  end

  Wx = diag(wx);
  W2x = sqrtm(Wx);
  W2xinv = inv(W2x);
  Ax = W2x*Cx*W2x;
  [Vx1,Dx1] = eig(Ax);
  Dx=diag(sort(diag(Dx1),'descend'));       % The following commands sort D and V from largest to smallest
  [c, ind]=sort(diag(Dx1),'descend');
  Vx = Vx1(:,ind);
%  [Vx,Dx] = sortem(Vx,Dx);                   % Sort D and V from largest to smallest
  Phix = W2xinv*Vx;

  Wy = diag(wy);
  W2y = sqrtm(Wy);
  W2yinv = inv(W2y);
  Ay = W2y*Cy*W2y;
  [Vy1,Dy1] = eig(Ay);
   Dy=diag(sort(diag(Dy1),'descend'));       % The following commands sort D and V from largest to smallest
  [c, ind]=sort(diag(Dy1),'descend');
  Vy = Vy1(:,ind);
  Phiy = W2yinv*Vy;

  W = diag(w);
  W2 = sqrtm(W);
  W2inv = inv(W2);
  A = W2*C*W2;
  [V2d,D2d] = eig(A);
  [V2d,D2d] = sortem(V2d,D2d);                   % Sort D and V from largest to smallest
  Phi2d = W2inv*V2d; 
  eigvals2d = diag(D2d);

  ind = 1;
  for i = 1:Nquadx
    for j = 1:Nquady
      D(ind,ind) = Dx(i,i)*Dy(j,j);
      ind = ind+1;
    end
  end
  clear ind

  Phi = zeros(Nquadx*Nquady,Nquadx*Nquady);
  for j = 1:Nquadx
    for i = 1:Nquadx
      Phixt(1+(i-1)*Nquady:i*Nquady,:) = Phiy*Phix(i,j);
    end
    Phi2(:,1+(j-1)*Nquady:j*Nquady) = Phixt;
  end

  [eigvals,ind] = sort(diag(D),'descend');
  Phi = Phi2(:,ind);

%
% Compute the convergence ratios
%

  for n = 1:N
    r(n) = sum(eigvals(1:n))/sum(eigvals(1:N));
  end

  figure(1)
  semilogy(eigvals/eigvals(1),'*','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('n')
  ylabel('\lambda_n/\lambda_1')

  figure(2)
  loglog(r,'*','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('n')
  ylabel('r_N')

%
% Compute the eigenfunctions and KL expansion
%

  KL = zeros(Nquadx,Nquady);
  for k = 1:Nmodes
    eigfun = zeros(Nquady,Nquadx);
    for j=1:Nquadx
      eigfun(:,j) = Phi(1+(j-1)*Nquady:j*Nquady,k);
    end
    Eigfun = eigfun';
    KL = KL + sqrt(eigvals(k))*Eigfun*randn(1);
  end

  KL2 = sqrt(eigvals(1))*Phix(:,1)*(Phiy(:,1)')*randn(1) ...
        + sqrt(eigvals(2))*Phix(:,1)*(Phiy(:,2)')*randn(1) ...
        + sqrt(eigvals(3))*Phix(:,2)*(Phiy(:,1)')*randn(1) ...
        + sqrt(eigvals(4))*Phix(:,2)*(Phiy(:,2)')*randn(1) ...
        + sqrt(eigvals(5))*Phix(:,1)*(Phiy(:,3)')*randn(1) ...
        + sqrt(eigvals(6))*Phix(:,3)*(Phiy(:,1)')*randn(1) ...
        + sqrt(eigvals(7))*Phix(:,1)*(Phiy(:,4)')*randn(1) ...
        + sqrt(eigvals(8))*Phix(:,4)*(Phiy(:,1)')*randn(1) ...
        + sqrt(eigvals(9))*Phix(:,2)*(Phiy(:,3)')*randn(1) ...
        + sqrt(eigvals(10))*Phix(:,3)*(Phiy(:,2)')*randn(1);

  for i = 1:Nquady
    X(:,i) = x;
  end
  for j = 1:Nquadx
    Y(j,:) = y';
  end

  figure(3)
  contourf(X,Y,KL)
  axis('equal')
  set(gca,'Fontsize',[20]);
  
