view +noname/calculateSolution.m @ 87:0a29a60e0b21

In Curve: Rearranged for speed. arc_length_fun is now a property of Curve. If it is not supplied, it is computed via the derivative and spline fitting. Switching to the arc length parameterization is much faster now. The new stuff can be tested with testArcLength.m (which should be deleted after that).
author Martin Almquist <martin.almquist@it.uu.se>
date Sun, 29 Nov 2015 22:23:09 +0100
parents 54d3ab296ba0
children 8298734b1938
line wrap: on
line source

% Calculates the solution of discretization for a given set of ms ts and orders.
%    discrHand -- function handle to a Discretization constructor
%    method    -- time stepping method
%    m         -- grid parameter
%    order     -- order of accuracy of the approximtion
%    T         -- time to calculate solution for
%    input paramters m, t, order may all be vectors.
function [] = calculateSolution(filename, discrHand, method, m, T_in, order, force_flag)
    default_arg('force_flag',false);

    if exist(filename,'file') && ~force_flag
        fprintf('File ''%s'' already exist.',filename);
        do_append = yesnoQuestion('Do you want to append to it?');
        if ~do_append
            fprintf('Exiting...\n');
            return
        end
    end

    sf = SolutionFile(filename);



    orderWidth = findFieldWidth('%d',order);
    mWidth = findFieldWidth('%d',m);
    TWidth = findFieldWidth('%d',T_in);

    for i = 1:length(order)
        for j = 1:length(m)
            T = sort(T_in); % Make sure times are sorted

            discr = discrHand(m(j),order(i));
            k_max = discr.getTimestep(method);

            % Do we want to to save the initial conditions?
            if T(1) == 0
                snapshot = discr.getTimeSnapshot(0);
                saveToFile(sf, method, order(i), m(j),T(1), snapshot, NaN, NaN, discr);
                T(1) = [];
            end

            % Find out if times to be calulated are integer multiples of the smallest one.
            time_multiples = T/T(1);

            is_int_multiples = all(time_multiples == int64(time_multiples));

            if is_int_multiples
                fprintf('Calculating time series in increments\n');
            else
                fprintf('Restarting for each time in timeseries\n');
            end

            % T now contains all the times we need to step to,
            % if T contained 0 it has now been removed.

            if is_int_multiples
                % Times are integer multiples, we can save time
                [k,N] = alignedTimestep(k_max,T(1));
                ts = discr.getTimestepper(method,k);
                runtime = 0;
                for l = 1:length(T)
                    end_step = N * time_multiples(l);
                    fprintf('[order = %-*d, m = %-*d, T = %-*d]: ',orderWidth,order(i),mWidth,m(j),TWidth,T(l));
                    clock_start = tic();
                    ts.stepN(end_step-ts.n,true);
                    runtime = runtime + toc(clock_start);
                    snapshot = discr.getTimeSnapshot(ts);
                    saveToFile(sf, method, order(i), m(j),T(l), snapshot, runtime, k, discr);
                    fprintf('Done! (%.3fs)\n',runtime);
                end
            else
                % Times are not interger multiples, we have to start from 0 every time.
                for l = 1:length(T)
                    [k,N] = alignedTimestep(k_max,T(l));
                    ts = discr.getTimestepper(method,k);
                    fprintf('[order = %-*d, m = %-*d, T = %-*d]: ',orderWidth,order(i),mWidth,m(j),TWidth,T(l));
                    clock_start = tic();
                    [v,t] = ts.stepN(N-ts.n,true);
                    runtime = toc(clock_start);
                    snapshot = discr.getTimeSnapshot(ts);
                    saveToFile(sf, method, order(i), m(j),T(l), snapshot, runtime, k, discr);
                    fprintf('Done! (%.3fs)\n',runtime);
                end

            end
            sf.stupidSave();
        end
    end
end


function saveToFile(sf, method, order, m, T, snapshot, runtime, k, discr)
    key.method = method;
    key.order  = order;
    key.m      = m;
    key.T      = T;

    entry.repr = snapshot;
    entry.runtime = runtime;
    entry.k = k;
    entry.discr = discr;

    sf.store(key,entry);
end