% Written by Vesa Kaarnioja, 2022
% Demonstration on how to solve parametric PDE source and eigenvalue
% problems over a Monte Carlo point set

% Precompute the FEM data
level = 5;
[mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(level);

% Set up the random diffusion coefficient
s = 100; % stochastic dimension
decay = 2; % the "decay parameter" \vartheta from the slides
indices = (1:s)';
deterministic = indices.^(-decay) .* sin(pi*indices*centers(:,1)') .* sin(pi* indices * centers(:,2)');
a = @(y) 2+deterministic'*y;

% Initialize the source term and loading vector
f = @(x) x(:,1); % Source term, in this case f(x) = x_1
rhs = mass(interior,:)*f(nodes); % Loading term; not affected by random diffusion coefficient

% Initialize the parallel pool (uses the Parallel Computing Toolbox)
poolobj = gcp('nocreate');
if isempty(poolobj)
    parpool;
end

null = zeros(ncoord,1);

% Part 1: Source problem
%% Sanity check: Can we solve the PDE for one realization of the uniform
% random vector y?
disp('Solve one PDE source problem')
mcnode = rand(s,1)-.5;
A = UpdateStiffness(grad,a(mcnode));
sol = null;
sol(interior) = A(interior,interior)\rhs;

% Visualize the results
figure;
trisurf(element,nodes(:,1),nodes(:,2),sol);

%% Solve PDE for many realizations of a random vector
n = 2^10;
disp('Solve many PDE source problems');
results = zeros(ncoord,n);
parfor ii = 1:n
    mcnode = rand(s,1)-.5;
    A = UpdateStiffness(grad,a(mcnode));
    sol = null;
    sol(interior) = A(interior,interior)\rhs;
    results(:,ii) = sol;    
end

% Visualize the results
for ii = [1,60,1024]
    figure;
    trisurf(element,nodes(:,1),nodes(:,2),results(:,ii));
end

% Part 2: Eigenvalue problem
%% Sanity check: Can we solve the eigenpair for one realization of the uniform
% random vector y?
disp('Solve one PDE eigenvalue problem');
mcnode = rand(s,1)-.5;
stiffness = UpdateStiffness(grad,a(mcnode));
[evecs,evals]= eigs(stiffness(interior,interior),mass(interior,interior),1,'smallestabs');
coef = sqrt(evecs' * mass(interior,interior) * evecs);
evecs = sign(evecs(1))*evecs/coef;
evecs_full = null;
evecs_full(interior) = evecs;

% Visualize the results
figure;
trisurf(element,nodes(:,1),nodes(:,2),evecs_full);

%% Solve the PDE for many realizations of the uniform random vector
n = 2^10;
disp('Solve many PDE eigenvalue problems');
eigenvalues = zeros(1,n);
eigenvectors = zeros(ncoord,n);
parfor ii = 1:n
    mcnode = rand(s,1)-.5;
    stiffness = UpdateStiffness(grad,a(mcnode));
    [evecs,evals] = eigs(stiffness(interior,interior),mass(interior,interior),1,'smallestabs');
    coef = sqrt(evecs' * mass(interior,interior) * evecs);
    evecs = sign(evecs(1))*evecs/coef;
    evecs_full = null;
    evecs_full(interior) = evecs;
    eigenvalues(ii) = evals;
    eigenvectors(:,ii) = evecs_full;
end

% Visualize the results
for ii = [1,60,1024]
    disp(['Minimum eigenvalue: ',num2str(eigenvalues(ii))]);
    figure;
    trisurf(element,nodes(:,1),nodes(:,2),eigenvectors(:,ii));
end
