Create residual viscosity package +rv and generalize the ResidualViscosity class - Generalize residual viscosity, by passing user-defined flux and calculating the time derivative outside of the update. - Create separate RungekuttaRV specifically using interior RV updates - Separate the artifical dissipation operator from the scheme AdvectionRV1D so that the same scheme can be reused for creating the diff op used by the ResidualViscosity class

% calculate a N-D mononomial with powers k in points x:
%  z = x(:,1).^k(1) * x(:,2).^k(2) * ...
function z = mononomial(x, k)
    assert(size(x,2) == length(k), 'k must have the same length as the width of x');

    if any(k < 0)
        z = x(:,1)*0;
        return
    end

    denom = prod(factorial(k));

    for i = 1:length(k)
        x(:,i) = x(:,i).^k(i);
    end
    z = prod(x,2)/denom;
end