%
% This tests the use of multilevel Monte Carlo for 
% Burgers equation with random initial data:
%
% PDE: du/dt + d/dx(1/2 u^2) = 0
%
% u(0,t) = 1, u(1,t)=-1,  t>=0
% u(x,0) =  U_1, 0<x<=0.5
% u(x,0) = -U_2, 0.5<x<1
% U_1, U_2 are independent uniform [0,1] random variables
%
% the output is P = 10 \int_0^1 u(x,0.2) dx
%

function hyp

close all; clear all;

%
% use global variables to pass parameters into gbm_l function
%

% global 

%
% call mlmc_test to perform MLMC tests
%

nvert = 3;      % plotting option (1 for slides, 2 for papers, 3 for full set)
M     = 4;      % refinement cost factor (2^gamma in general MLMC Thm)

N     = 1000;   % samples for convergence tests
L     = 8;      % levels for convergence tests 

N0    = 100;   % initial number of samples on each MLMC level
Eps   = [ 0.005 0.01 0.02 0.05 0.1 ];  % desired accuracies for MLMC calcs

figs = mlmc_test(@hyp_l, M, N,L, N0,Eps, nvert);


%-------------------------------------------------------
%
% level l estimator for Burgers equation solver
%

function [sum1 sum2] = hyp_l(l,N)

% global

r  = 0.25;   % timestep/spacing
c  = 0.1;    % numerical smoothing coefficient

nf = 2^(l+1);
hf = 1/nf;

if l>0
  nc = nf/2;
  hc = 1/nc;
end

sum1(1:4) = 0;
sum2(1:2) = 0;

for N1 = 1:100:N
  N2 = min(100,N-N1+1);
  U  = rand(2,N2);

  uf = [ ones(1,N2); ...         % x=0
         ones(nf/2,1)*U(1,:);    % 0<x<=0.5
        -ones(nf/2-1,1)*U(2,:);  % 0.5<x<1
        -ones(1,N2) ];           % x=1

  i  = 2:nf;

  for n = 1:nf                         % time-marching
    uf(i,:) = uf(i,:) - 0.25*r*(uf(i+1,:).^2-uf(i-1,:).^2) ...
                      +      c*(uf(i+1,:)-2*uf(i,:)+uf(i-1,:));
  end

  Pf = 10*sum(hf*uf.^2,1);                % output

  if l==0
    Pc = zeros(1,N2);
  else

    uc = [ ones(1,N2); ...         % x=0
           ones(nc/2,1)*U(1,:);    % 0<x<=0.5
          -ones(nc/2-1,1)*U(2,:);  % 0.5<x<1
          -ones(1,N2) ];           % x=1

    i  = 2:nc;

    for n = 1:nc                       % time-marching
      uc(i,:) = uc(i,:) - 0.25*r*(uc(i+1,:).^2-uc(i-1,:).^2) ...
                        +      c*(uc(i+1,:)-2*uc(i,:)+uc(i-1,:));
    end

    Pc = 10*sum(hc*uc.^2,1);              % output
  end

  sum1(1) = sum1(1) + sum(Pf-Pc);
  sum1(2) = sum1(2) + sum((Pf-Pc).^2);
  sum1(3) = sum1(3) + sum((Pf-Pc).^3);
  sum1(4) = sum1(4) + sum((Pf-Pc).^4);
  sum2(1) = sum2(1) + sum(Pf);
  sum2(2) = sum2(2) + sum(Pf.^2);
end
