%            Aysajan Abidin

% Function:  U = Solve_spde(u0,t0,T,K,H,S,NOISE)
%            Solves parabolic spdes driven by additive noise on rectangular
%            domains in 1D for different spatial discretization.
%
% Input  
%      H:    The number of space steps in the finest spatial partition. So
%            we have H-1 points inside the interval (0,1).
%      K:    The number of time steps in the finest time partition.
%            The length of each time step is dt = T/K, where T is 
%            described below. 
%            The finest means the best possible space and time partitions
%            for which the specific computer can affort the computational
%            power. In our test we take H = 2^11 and K = 2^11
%      T:    Final timepoint.
%     u0:    Initial Contidion
%     t0:    Initial timepoint.
%      S:    A vector used to specifie the length of each spatial interval.
%            Example: S = [1,4,8,16,32,64,128] means that there are 7
%            spatial intervals of different length, each of length
%            Dh=S(i)*dh, where dh=1/H. Note that mod(H,S(i))=0, for all i.
%  NOISE:    Specifies whether the noise is White or Colored
% 
% Output
%        U:  A matrix consisting in each column of the solution to the SDPE at the final
%            time step for each spatial discretization.
 
    function U = Solver_spde(u0,t0,T,K,H,S,NOISE)
    
    %%%%%%% look for errors in input:
    if mod(H,S) ~= 0
        error('Error! Note that mod(H,S(i))=0 for all i.');
    end
    if and(strcmp(NOISE,'White')==0,strcmp(NOISE,'Colored') == 0)
        error('Error! Type of noise wrongly specified');
    end
    
    %%%%%%% Setting parameters related to the time and the space domain
    hh = length(S);          % Number of different spatial partitions
    L = H./S;                % Number of space steps in each spatial partition
    dh = 1/H;                % Length of space step in the finest spatial partition.
    dt = (T-t0)/K;           % Length of time step in the finest time partition.
    %%%%%%% Generation of Brownian paths
    % dB is a 2-d Matrix of 1D Brownian increments (delta beta). 
    % A different Brownian path in each line, different time point in 
    % each column.  Note that the Brownian increments are not space dependent.
    % These Brownian increments relate to the basic time step, dt.
    dB = sqrt(dt)*randn(H-1,K);
   
    B = cumsum(dB,2);    % B is a matrix of the Brownian paths corresponding to dB
    dB_step = diff([zeros(H-1,1) B(:,1:end)],1,2);  %Calculates the Brownian increment
    U=[]; W_org=[]; 
    for i = 1:hh % Loop over different spatial partitions
        % Create the parametes for the assembly of stiffness and mass matrix
        p = 1/L(i):1/L(i):(1-1/L(i));
        n = length(p);
        t = [1:(n-1); 2:n; ones(1,(n-1))];
        e = [1, n; 1,  2];
       [stiff mass] = stiff_mass(p,e,t);  % Assemble stiffness and mass matrices  
       disp(sprintf('subround %d', i))       % printing i to know how the process is progressing.
       message = sprintf('  Mass & Stiffness Matrices are ready.');
            disp(message);
       
       DB = dB_step(S(i):S(i):end,:); % Take the corresponding Brownian paths for different space steps
       
       if strcmp(NOISE, 'White')
           Lower = chol(mass)';           
       else if strcmp(NOISE, 'Colored')               
               Lower = Assemble_noise(p); %Assemble the cholesky factor.
           end                     
       end
       message = sprintf('  Cholesky factor is ready. \n');
               disp(message);       
       W_org = Lower * DB;           % Assemble the noise

       u_prev = u0(S(i):S(i):end);    
       for jj = 1:K              
           W = W_org(:,jj);     
           % Compute the finite element solution at time t = dt
           u_prev = (mass+stiff*dt)\(mass*u_prev+W);                    
       end
       if length(u_prev) == H-1
           U = [U u_prev];
       else
           N = H - length(u_prev);
           U(:,i) = [u_prev; zeros(N-1,1)];
       end
    end