% Driver for the one-dimensional labs.
%
% The editable script supplies two functions:
%
%     names = method_names()                 a cell array of method labels
%     [x, u] = solve1d(mid, N, f, ua, ub)    method mid, resolution N
%
% solve1d picks its own grid x (uniform, graded, Chebyshev, whatever the method
% wants), returns the solution values on that grid, and the driver measures the
% error against the manufactured solution at exactly those points. f is a
% function handle giving the source term, so a method is free to evaluate it
% anywhere.
%
% Timings come from tic/toc around the solve and nothing else. Each method gets
% a warm-up call before its first timed one, because the first call also pays
% for compiling the method, and anything under 5 ms is repeated with the best
% kept: a sub-millisecond solve otherwise measures the clock.
p = lp_params();
kind = p.kind;
names = method_names();
nm = numel(names);
ub = lp_exact1d(kind, [0; 1]);
f = @(xx) drv_rhs(kind, xx);

out = struct();
out.names = names;
out.kind = kind;

if strcmp(p.mode, 'single')
    t0 = tic;
    [x, u] = solve1d(p.mid, p.N, f, ub(1), ub(2));
    out.ms = toc(t0) * 1000;
    uex = lp_exact1d(kind, x(:));
    err = u(:) - uex(:);
    out.x = x(:)';
    out.u = u(:)';
    out.uex = uex(:)';
    out.err = err';
    out.maxerr = max(abs(err));
    out.n = numel(x) - 1;
    xf = linspace(0, 1, p.ngrid)';
    out.xf = xf';
    out.uf = lp_exact1d(kind, xf)';
else
    ns = p.ns(:)';
    nn = numel(ns);
    E = NaN(nm, nn);
    T = NaN(nm, nn);
    Q = NaN(nm, nn);
    nmax = p.nmaxs(:)';
    timed = p.timed;
    for k = 1:nm
        [xw, uw] = solve1d(k, 8, f, ub(1), ub(2));   % warm-up
        for j = 1:nn
            N = ns(j);
            if N > nmax(k)
                continue;
            end
            t0 = tic;
            [x, u] = solve1d(k, N, f, ub(1), ub(2));
            el = toc(t0);
            rep = 0;
            while timed > 0 && el < 0.005 && rep < 4
                rep = rep + 1;
                t0 = tic;
                [x, u] = solve1d(k, N, f, ub(1), ub(2));
                el = min(el, toc(t0));
            end
            uex = lp_exact1d(kind, x(:));
            E(k, j) = max(abs(u(:) - uex(:)));
            T(k, j) = el * 1000;
        end
        Q(k, :) = lp_orders(E(k, :), ns);
    end
    out.ns = ns;
    out.E = lp_rows(E);
    out.T = lp_rows(T);
    out.orders = lp_rows(Q);
end

lp_write(out);

function y = drv_rhs(kind, x)
[~, y] = lp_exact1d(kind, x);
end