% Written by Vesa Kaarnioja, 2023
% Demonstration on how to compute a 2048-dimensional reference solution, to be used for dimension truncation error estimation

clear all; close all;

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

% Set up the diffusion coefficient
s = 2048; % stochastic dimension
decay = 2.0; % decay of the input random field
indices = (1:s)';
deterministic = indices.^(-decay) .* sin(pi*indices*centers(:,1)') .* sin(pi* indices * centers(:,2)');
a = @(y) 2 + deterministic'*y;

% Set up the loading
f = @(x) x(:,1);
rhs = mass(interior,:)*f(nodes);

% Import the generating vector
gen = load('offtheshelf2048.txt');
z = gen(1:s); % truncate generating vector to s terms
n = 2^15; % fix number of QMC nodes
one = ones(length(interior),1);

% Solve 2^15 PDEs in parallel
tmp = 0;
parfor ii = 1:n
    qmcnode = mod(z*ii/n,1)-1/2;
    A = UpdateStiffness(grad,a(qmcnode));
    sol = A(interior,interior)\rhs;
    tmp = tmp + one'*mass(interior,interior)*sol;
end

% Compute the QMC average
reference = tmp/n;
disp(reference)
