function dgclaw_compute(par)
%DGCLAW_COMPUTE
%   Discontinuous Galerkin solver for systems of hyperbolic
%   conservation laws in 1D.
%
%   dimension 1: system component
%   dimension 2: x
%   dimension 3: polynomial moments
%
%   2012/06/18 by Benjamin Seibold

%================================================================
% Parameters
%================================================================
try par.name;          catch, par.name = ''; end
try par.len;           catch, par.len = 1; end
try par.nx;            catch, par.nx = 100; end
try par.np;            catch, par.np = 2; end
try par.cfl;           catch, par.cfl = .99/(2*par.np+1); end
try par.M;             catch, par.M = 0; end
try par.uax;           catch, par.uax = [-1 1]*1e-6; end
try par.evolve;        catch, par.evolve = nan; end
try par.g;             catch, par.g = nan; end
try par.bc;            catch, par.bc = nan; end
try par.tfinal;        catch, par.tfinal = 1; end
try par.steps;         catch, par.steps = 100; end
try par.tplot; par.steps = length(par.tplot)-1;
catch, par.tplot = linspace(0,par.tfinal,par.steps+1); end
try par.flag_plot;     catch, par.flag_plot = 0; end
try par.flag_savefig;  catch, par.flag_savefig = 0; end
try par.flag_savefile; catch, par.flag_savefile = 0; end
try par.savepath;      catch, par.savepath = '.'; end

%================================================================
% Initialization
%================================================================
par.h = par.len/par.nx;
x = par.h/2:par.h:par.len-par.h/2;
par.maxcurv = par.M*par.h^2; % for relaxed limiting
U = limit(par,moments_ic(par,x),0); % initial conditions
n_add = size(U,1)-size(par.uax,1);
par.uax(end+1:end+n_add,:) = ones(n_add,1)*[-1 1]*1e-6;

%================================================================
% Time loop
%================================================================
t = 0;
for j = 0:par.steps % loop over output steps
    par.plot_step = j;
    t_new = par.tplot(par.plot_step+1); % new time to reach
    U = evolve(par,x,U,[t,t_new]);
    t = t_new;
    par = visualize(par,x,U,t); % visualize new state
end

%================================================================
% Computational functions
%================================================================
function Unew = evolve(par,x,U,tstep)
% Evolve solution from tstep(1) to tstep(2) with as many
% computation steps it takes.
while diff(tstep)>0
    lambda = max(max(abs(...    % maximum characteristic
        par.feigs(par,U(:,:,1),x,tstep(1)))));  % velocity
    dtmax = par.cfl*par.h/lambda; % admissible time step
    dt = min(dtmax,diff(tstep)); % actually used time step
    U = shu_osher(par,U,tstep(1)+[0 dt]); % perform ODE step
    tstep(1) = tstep(1)+dt; % update current time
end
Unew = U;

function par_new = visualize(par,x,U,t)
x = [x(1)-par.h,x,x(end)+par.h];
U = add_ghost_cells(par,U,t);
[px,pu] = numerical_solution(x,U);
ns = size(U,1);
sax = ceil(sqrt(ns)); sax = [sax,ceil(ns/sax)];
if (par.flag_savefig||par.flag_savefile)&&~exist(par.savepath,'dir')
    mkdir(par.savepath) % create sub-directory if needed
end
if par.flag_plot||par.flag_savefig % plot visibly or save figure
    if ~par.flag_plot
        figure('visible','off') % open invisible figure window
    end
    clf
    for js = 1:ns % loop over system components
        par.uax(js,1) = min(par.uax(js,1),min(U(js,:,1),[],2));
        par.uax(js,2) = max(par.uax(js,2),max(U(js,:,1),[],2));
        subplot(sax(1),sax(2),js)
        plot(px,pu(js,:),'r-','linewidth',1.5)
        axis([0 par.len par.uax(js,:)+diff(par.uax(js,:))*[-1 1]*.1])
        title(sprintf('%s: u[%d] at t=%0.2f',par.name,js,t))
        pos = get(gca,'position');
        set(gca,'position',[.05 pos(2) .93 pos(4)])
    end
    if par.flag_savefig
        % save figure into png file
        fname = sprintf('%s/fig_%06d.png',par.savepath,par.plot_step);
        res = 800;
        set(gcf,'paperpositionmode','auto')
        set(gcf,'position',[10 100 res res*.8])
        print(gcf,'-dpng','-r290',fname);
    else
        drawnow % plot without saving
    end
    if ~par.flag_plot
        close(gcf) % close invisible figure window
    end
end
if par.flag_savefile % save solution into file
    fname = sprintf('%s/savefile_%06d',par.savepath,par.plot_step);
    save(fname,'par','x','U','t')
end
par.plot_step = par.plot_step+1;
par_new = par;

function U = moments_ic(par,x)
% Compute moments of initial conditions with respect to Legendre
% polynomials up to order par.np.
x = x(:); % make x a column vector
nx = length(x);
h = x(2)-x(1);
% use fourth dimension for Gauss points
[xg,wg] = gauss_points; ng = length(xg);
xgh = h/2*xg'; % cell-scaled Gauss points
X = permute(x(:,ones(ng,1))+permute(xgh(:,ones(nx,1)),[2 1]),...
    [3 1 4 2]); % in-cell-Gauss-point grid
f = par.ic(par,X);
F = f(:,:,ones(par.np+1,1),:); % [ns,nx,np+1,ng]
ns = size(F,1);
hi = (2*(0:par.np)'+1)/2; % half integers
P = leg_poly(xg,0:par.np).*hi(:,ones(ng,1));
P = permute(P(:,:,ones(ns,1),ones(nx,1)),[3 4 1 2]);
U = wsum(F.*P,wg,4);

function Unew = shu_osher(par,U,tstep)
% Shu-Osher 3rd order SSP step from tstep(1) to tstep(2),
% appropriately limited.
k = diff(tstep);
U1 = limit(par,euler(par,U,tstep(1),k),tstep(2));
U2 = limit(par,3/4*U+1/4*euler(par,U1,tstep(2),k),tstep(2));
Unew = limit(par,1/3*U+2/3*euler(par,U2,tstep(1)+k/2,k),tstep(2));

function Unew = euler(par,U,t,k)
% Forward Euler step with U of size k, with rhs evaluated at t
Unew = U+k*rhs(par,U,t); % forward Euler step
if isa(par.evolve,'function_handle')
    Unew = par.evolve(Unew,t,k); % additional evolution
end

function dU = rhs(par,U,t)
% Right hand side of state U, at time t.
[ns,nx,np] = size(U);
x = par.h/2:par.h:par.len-par.h/2;
% fluxes through left and right cell boundaries
U_aug = add_ghost_cells(par,U,t);
ul = wsum(U_aug,(-1).^(0:np-1),3); % values of left side of cells
ur = wsum(U_aug,ones(1,np),3);     % values on right side of cells
fluxl = f_num(par,ur(:,1:end-2),ul(:,2:end-1),x-par.h/2,t);
fluxr = f_num(par,ur(:,2:end-1),ul(:,3:end)  ,x+par.h/2,t);
FL = fluxl(:,:,ones(np,1));
FL(:,:,2:2:end) = -FL(:,:,2:2:end);
FR = fluxr(:,:,ones(np,1));
% integrals over cells
Q = quadrature_integrals(par,U,t);
% putting everything together
oddind = 2*(0:np-1)'+1; % odd indices
C = permute(oddind(:,ones(ns,1),ones(nx,1)),[2 3 1]);
dU = C.*(Q-(FR-FL))/par.h;

function Y = quadrature_integrals(par,U,t)
% Gauss-Legendre quadrature of integral terms at time t.
% This is all local, without communication between cells.
[ns,nx,np] = size(U);
x = (par.h/2:par.h:par.len-par.h/2)';
% use fourth dimension for Gauss points
[xg,wg] = gauss_points; ng = length(xg);
xgh = par.h/2*xg'; % cell-scaled Gauss points
X = permute(x(:,ones(ng,1)),[3 1 2])+...
    permute(xgh(:,ones(nx,1)),[3 2 1]); % positions
p = leg_poly(xg,0:np-1); % Legendre polynomials
P = permute(p(:,:,ones(ns,1),ones(nx,1)),[3 4 1 2]); % [ns,nx,np,ng]
p1 = leg_poly(xg,0:np-1,1); % derivatives of Legendre polynomials
P1 = permute(p1(:,:,ones(ns,1),ones(nx,1)),[3 4 1 2]); % [ns,nx,np,ng]
u = permute(reshape(p'*reshape(permute(U,[3 2 1]),np,[],1),...
    ng,nx,[]),[3 2 1]); % u-values at Gauss points
fu = par.f(par,u,X,t); % evaluate hyperbolic flux function
FU = permute(fu(:,:,:,ones(np,1)),[1 2 4 3]); % [ns,nx,np,ng]
gu = par.g(par,u,X,t); % evaluate source term of balance law system
GU = permute(gu(:,:,:,ones(np,1)),[1 2 4 3]); % [ns,nx,np,ng]
Y = wsum(FU.*P1+par.h/2*GU.*P,wg,4);

function z = f_num(par,a,b,x,t)
% Numerical flux function of states a and b, evaluated
% at position x and time t.
% Here: local Lax-Friedrichs
la = par.feigs(par,a,x,t);
lb = par.feigs(par,b,x,t);
ns = size(la,1);
m = max(abs([la;lb])); % maximal characteristic velocity
c = m(ones(ns,1),:); % [ns,nx]
fa = par.f(par,a,x,t);
fb = par.f(par,b,x,t);
z = (fa+fb-c.*(b-a))/2;

function Unew = limit(par,U,t)
% Apply limiters to reconstructed solution.
% Time matters for boundary conditions.
[ns,nx,np] = size(U);
Unew = U;
if np<=1, return, end % no limiting required
sl = wsum(U,[0,(-1).^(0:np-2)],3); % slope at left
sr = wsum(U,[0,ones(1,np-1)],3); % slope at right
sl_lim = limited_slope(par,U(:,:,1),sl,t);
sr_lim = limited_slope(par,U(:,:,1),sr,t);
s_lim = limited_slope(par,U(:,:,1),U(:,:,2),t);
%U_lim = cat(3,U(:,:,1),s_lim,zeros(ns,nx,np-2)); % more restrictive
U_lim = cat(3,U(:,:,1),s_lim/2,zeros(ns,nx,np-2)); % less restrictive
ind = (abs(sr_lim-sr)+abs(sl_lim-sl))>1e-8; % different result
ind = ind(:,:,ones(np,1));
Unew(ind) = U_lim(ind);

function v = limited_slope(par,u0,u1,t)
% Limit slopes given in u1 vs. data given in u0.
% Both inputs must be of length 1 in dimension 3.
% Time matters for boundary conditions.
u_aug = add_ghost_cells(par,u0,t);
u0_l = u0-u_aug(:,1:end-2); % difference to left
u0_r = u_aug(:,3:end)-u0;   % difference to right
w1 = par.fL(par,u0,u1); % transform cell slope into char. variables
w0_l = par.fL(par,u0,u0_l); % transform left-side slope into char. var.
w0_r = par.fL(par,u0,u0_r); % transform right-side slope into char. var.
w = minmod_tvb(cat(3,w1,w0_l,w0_r),par.maxcurv);
v = par.fR(par,u0,w); % transform limited slope back

function Unew = add_ghost_cells(par,U,t)
% Augments solution array by ghost cells left and right of
% computational domain.
if isa(par.bc,'function_handle')
    [ns,~,np] = size(U);
    e = zeros(ns,1,np-1);
    [u0,u1] = par.bc(t);
    Unew = [cat(3,u0,e),U,cat(3,u1,e)]; % Dirichlet b.c.
else
    Unew = [U(:,end,:),U,U(:,1,:)]; % periodic b.c.
end

%================================================================
% Auxilliary functions
%================================================================
function [px,pu] = numerical_solution(x,U)
% Generate numerical solution function for plotting.
nc = 8; % number of evaluation points in each cell
[ns,nx,np] = size(U);
h = x(2)-x(1);
xc = linspace(-1,1,nc);
p = leg_poly(xc,0:np-1);
px = reshape(ones(nc,1)*x+h/2*xc'*ones(1,nx),1,[]);
pu0 = p'*reshape(permute(U,[3 2 1]),np,[],1);
pu = permute(reshape(pu0,1,[],ns),[3 2 1]);

function y = leg_poly(x,deg,ndiff)
% Evaluates Legendre polynomials of orders encoded in
% vector deg. If ndiff is provided, then ndiff-th derivative
% is taken.
if nargin<3, ndiff = 0; end
x2 = x.^2;
n = length(deg);
y = zeros(n,length(x));
if ndiff==0
    for j = 1:n
        switch deg(j)
            case 0, y(j,:) = ones(size(x));
            case 1, y(j,:) = x;
            case 2, y(j,:) = (3*x2-1)/2;
            case 3, y(j,:) = (5*x2-3).*x/2;
            case 4, y(j,:) = ((35*x2-30).*x2+3)/8;
            case 5, y(j,:) = ((63*x2-70).*x2+15)*x/8;
        end
    end
elseif ndiff==1
    for j = 1:n
        switch deg(j)
            case 0, y(j,:) = zeros(size(x));
            case 1, y(j,:) = ones(size(x));
            case 2, y(j,:) = 3*x;
            case 3, y(j,:) = (15*x2-3)/2;
            case 4, y(j,:) = (35*x2-15).*x/2;
            case 5, y(j,:) = ((315*x2-210).*x2+15)/8;
        end
    end
end

function [x,w] = gauss_points()
% Gauss-Legendre quadrature points and weights.
x1 = sqrt(5-2*sqrt(10/7))/3;
x2 = sqrt(5+2*sqrt(10/7))/3;
w1 = (322+13*sqrt(70))/900;
w2 = (322-13*sqrt(70))/900;
x = [-x2 -x1 0 x1 x2];
w = [w2 w1 128/225 w1 w2];

function v = minmod_tvb(U,maxcurv)
% Generalized minmod limiter.
v = minmod(U);
u0 = U(:,:,1);
ind = abs(u0)<=maxcurv;
v(ind) = u0(ind);

function v = minmod(U)
% Apply minmod operation over third dimension of U.
S = sign(U);
v = (S(:,:,1).*min(abs(U),[],3)).*(all(S>0,3)|all(S<0,3));

function B = wsum(A,w,dim)
% Dot-multiply array with vector w along dimension dim.
% This is essentially sum(A,dim), weighted with w.
if length(w)==1, B = w*A; return, end
dims = size(A);
ind = 1:length(dims); ind([1 dim]) = [dim 1];
sdims = dims(ind);
S = w(:)'*reshape(permute(A,ind),sdims(1),[]);
B = ipermute(reshape(S,[1 sdims(2:end)]),ind);
