function ts = getRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    switch opt.method
        case 'rkRv'
            ts = rkRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvBdf'
            ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvMs'
            ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvMg'
            ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvInstage'
            ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        otherwise
            error('Timestepping method ''%s'' not supported',method);
    end
end

function F = stabilizedRhs(D, data)
    if isa(data, 'function_handle')
        F = @(v, t, viscosity) D(v, viscosity) + data(v);
    else
        F = @(v, t, viscosity) D(v, viscosity) + data;
    end
end

function F = unstabilizedRhs(D, data)
    if isa(data, 'function_handle')
        F = @(v, t) D(v) + data(v);
    else
        F = @(v, t) D(v) + data;
    end
end

function ts = rkRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    ts = rv.time.RungekuttaRv(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end

function ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    ts = rv.time.RungekuttaRvBdf(F, opt.k, t0, v0, residualViscosity, opt.rkOrder, opt.bdfOrders);
end

function ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    F_unstab = unstabilizedRhs(diffOpStruct.D_unstable, data);
    ts = rv.time.RungekuttaRvMultiStage(F, F_unstab, opt.k, t0, v0,...
                                        residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end

function ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    F_coarse = unstabilizedRhs(diffOpStruct.D_coarse, data);
    ts = rv.time.RungekuttaRvMultiGrid(F, F_coarse, opt.k, t0, v0,...
                                       residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end

function ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    ts = rv.time.RungekuttaRvInstage(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end