%                 May 21st, 2007 -- Aysajan Abidin
%
%                 January 11th 2006 - Asta Herdis Hall
% Function:       [U U_T] = Solver_parabolic_spde(u0,t0,T,K,R,J,p,e,t,TAG)
%                 Solves parabolic spdes on rectangular domains in 1D or 2D
%                 depending on whether TAG=1 or 2 for different time partitions.
%
% Uses:           User defined function ei_m_l = Eigen_rec(lth_x,lth_y,Np);
%                 The code for the 2D case should be modified based on the
%                 new technique for assembling the noise presented in the
%                 thesis.
%              
% Input  
%           u0:   Initial conditions, i.e. an approximation of the 
%                 solution to the pde at time t=0.
%            J:   Indicates how many terms we use in the sum for
%                 approximating the Wiener process. It equals the number of
%                 spatial mesh points.
%        NOISE:   Takes values 'add' or 'mult', indicating whether
%                 we have an additive or multipicative noise term. Note
%                 that in 1-D it is always 'add', but that parameter has no
%                 effect in 1-D case of the code, since we direcly wrote
%                 the assembly code for the noise in the program.
%     alfa,u0,   
% K,T,t0,R,TAG:   See parameters_description.m
%        p,e,t:   See parameters_description.m
% Output
%        U:       The solution to the spde for all timepoints.
%        U_T:     The solution to the spde at time T.
 
function [U, U_T] = Solver_parabolic_spde(alfa,u0,t0,T,K,R,J,p,e,t,NOISE,TAG)

    %%%%%%% looking for errors in input:
    if max(abs(rem(K,R))) ~= 0
        error('Error in the specification of R');
    end
    if and(strcmp(NOISE,'mult')==0,strcmp(NOISE,'add') == 0)
        error('Error - Type of noise wrongly specified');
    end
    
    %%%%%%% Setting parameters related to the time and the space domain
    kk = length(R);          % Number of different time partitions that we
                             % compute a solution for.
    L = K./R;                % Number of time steps in each time partition
    dt = (T-t0)/K;           % Length of time step in finest time partition.

    %%%%%%% Setting other parameters
    Np = size(p,2);          % Number of points in space mesh

    %%%%%%% 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(J,K);
    % B is a matrix of the Brownian paths corresponding to dB
    B = cumsum(dB,2);

if TAG==1    
    %%%%%%% Assemble stiffness- and mass matrices and boundary conditions 
    [stiff mass] = Assemble_stiff_mass(p,e,t,TAG);   
    message = sprintf('  Mass & Stiffness Matrices are ready.');
            disp(message);
    %%%%%%%
    U=[];U_T=[];W_org=[];
    for k = 1:kk       % Loop over different time steps
        U = [U u0];
        u_prev = u0;
        Dt = dt*R(k);   % Length of current time step, 
                        % a multiple of basic time step                        
        
        disp(sprintf('subround %d', i))       % printing i to know how the process is progressing.                
        
        % Find Brownian increment corresponding to current division
        % of the time interval    
        dB_step = diff([zeros(Np,1) B(:,R(k):R(k):end)],1,2);      
        
        Lower = chol(mass)';
        message = sprintf('  Cholesky factor is ready. \n');
               disp(message);  
        
        W_org = Lower*dB_step;
        
        for jj = 1:L(k)  % Time stepping, L(k) is the number of intervals
                         % in the current division of the time interval
            % Use backward Euler method to find next approximation of U
            
            W = W_org(:,jj);

            % Compute the finite element solution at time t = jj*Dt
            u_prev = (mass+stiff*Dt)\(mass*u_prev+W);

            U = [U u_prev];
        end
        U_T = [U_T u_prev];        
    end

else
    lth = max(p,[],2)-min(p,[],2);    % Length of x- and y-axis of domain
    %%%%%%% Computation of eigenvalues and eigenfunctions
    % Start by computing the eigenvalues of the Laplacian ...
    ei_m_l = Eigen_rec(lth(1),lth(2),J+1);
    % Do not use the first eigenvalue, which is 0
    ei_m_l = ei_m_l(2:end,:);
    % ... and the corresponding eigenvalues of Q, the covariance
    % matrix for the Wiener process:
    gamma_j = ei_m_l(:,1).^(-alfa);

    % Find eigenfunctions e_j, j = 1:J, in each of the mesh points:
    % Note that we have different eigenfunctions in each line and
    % different points in space in each column.
    tmp1 = sin(pi*repmat(ei_m_l(:,2),1,Np).*repmat(p(1,:)',1,J)'/lth(1));
    tmp2 = sin(pi*repmat(ei_m_l(:,3),1,Np).*repmat(p(2,:)',1,J)'/lth(2));
    e_n = 2/sqrt(lth(1)*lth(2))*(tmp1.*tmp2);

    %%%%%%% Assemble stiffness- and mass matrices and boundary conditions 
    [stiff mass] = Assemble_stiff_mass(p,e,t,TAG);   
                                  % One-time assembly of stiffness- a
                                  % and mass matrices.
    [Q,G] = Assemble_boundary(p,e,0); 
                                  % Boundary matrices.  Here, we assume
                                  % that the boundary conditions are
                                  % not time dependent.
                                
    %%%%%%%
    U=[];U_T=[];W_org=[];
    for k = 1:kk       % Loop over different time steps
        U = [U u0];
        u_prev = u0;
        Dt = dt*R(k);   % Length of current time step, 
                        % a multiple of basic time step                        
                        
        % Find Brownian increment corresponding to current division
        % of the time interval
    
        dB_step = diff([zeros(Np,1) B(:,R(k):R(k):end)],1,2);

        f_mat = Assemble_load(t0,T,K/R(k),p,e,t);  % Assemble time dep. 
                                                   % load vector                                                   

        % if we have additive noise then we assemble the W-vector here.
        % Note: we get a Np-x-M matrix where each column stands for a 
        % particular time point.
        if strcmp(NOISE,'add')
            W_org = Assemble_noise_additive(gamma_j,e_n,dB_step,t0,T,p,e,t);
        end
        for jj = 1:L(k)  % Time stepping, L(k) is the number of intervals
                         % in the current division of the time interval
            % Use backward Euler method to find next approximation of U
            t_prev = (jj-1)*Dt;
      
            % If we have multiplicative noise then we need to assemble the
            % W-vector in every time step.  
            if strcmp(NOISE,'mult')
                W = Assemble_noise_multiplicative(gamma_j,e_n,dB_step(:,jj),t_prev,u_prev,p,e,t);
            else
                W = W_org(:,jj);
            end
            %    [Q,G] = Boundary(p,e,jj*Dt); % uncomment if boundary matrices are
                                    % time-dependent.

            % Compute the finite element solution at time t = jj*Dt
            u_prev = (mass+(Q+stiff)*Dt)\(mass*u_prev+Dt*(G+f_mat(:,jj))+W);

            U = [U u_prev];
        end
        U_T = [U_T u_prev];
    end
end