%                 Aysajan Abidin

clf, clear all;
%%%%%%%% User specified parameters: %%%%%%%%%%%%
alfa = 1;
TAG = 1;
N = 100; % Number of simulations
K = 2^11;        
T = 1;
t0 = 0;
R = [1,2,4,8,16,];%,32,64,128,256];

% generate uniform mesh
p = 1/256:1/256:(1-1/256);%1/2048:1/2048:(1-1/2048);%fine space mesh
Np = length(p);
t = [1:(Np-1); 2:Np; ones(1,(Np-1))];
e = [1, Np;
     1, 2];

u0 = zeros(Np,1);

dt = (T-t0)/K; % Length of time step in finest time partition.
L = K./R;          % Number of time steps in each time partition

% Computing the finite element solutions and
% error estimates:
U=[]; U_T = []; error = [];
L2_err = 0;
for i = 1:N
    clear U, U_T;
    disp(sprintf('Round %d:', i))  % printing i to know how the process is progressing.
    [U U_T] = Solver_parabolic_spde(alfa,u0,t0,T,K,R,Np,p,e,t,'add',TAG);  
    L2_err = L2_err + L2_ErrorEst(U_T(:,1),U_T(:,2:end),p,e,t,TAG).^2;
end
error = [error ; sqrt(L2_err/N)] % computes L_2(Omega,H) error.

%% plot of the error on loglog scale
k = 1./L;
loglog(k(2:end),error);
hold on,
loglog(k(2:end),error,'r*');
xlabel('Length of the time step k'),ylabel('L_2(Omega,H) error for corresponding time steps')
save error;
