function diffOpsStruct = constructDiffOpsMultiGrid(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for stabilized solution vector
    [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    % DiffOp for unstabilized solution vector
    D_coarse = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    %% DiffOps for residual viscosity
    [D_t, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
    D_flux = @(v) -D_t(v);
    diffOpsStruct = struct('D_scheme', D_scheme,...
                     'D_coarse', D_coarse,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', penalties_scheme,...
                     'penalties_res', penalties_res);
end

% TODO: Not clear if this is needed
% TODO: Only works for equidistant grids
function [D_coarse, penalties_coarse] = constructCoarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
    m = g.m();
    lim = g.lim();
    m_coarse = (m-1)/2 + 1;
    g_coarse = grid.equidistant(m_coarse, lim{1});
    [D_coarse, penalties_coarse] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g_coarse, schemeOrder, schemeParams, opSet, BCs);
    x = g.x{1};
    x_coarse = x(1:2:end);
    % TODO: Fix interpolation
    D_coarse = @(v) interp1(x_coarse,D_coarse(v(1:2:end)),x,'spline');
end