% % 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]);