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

  clear all

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

  Lx = 0.5;
  Ly = 0.25;

  Nmodes = 104;
  Nquadx = 100;
  Nquady = 50;
  N = Nquadx*Nquady;
  hx = 2/Nquadx;
  hy = 1/Nquady;
  wx = hx*ones(1,Nquadx);
  wy = hy*ones(1,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

  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;

  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(1:1000)/eigvals(1),'*','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('n')
  ylabel('\lambda_n/\lambda_1')

  figure(2)
  loglog(r(1:1000),'*','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]);