function [D_scheme, D_coarse, D_flux, D_t, penalties_scheme, penalties_res] = constructDiffOpsMultiGrid(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
    % DiffOps for solution vector
    [D_scheme, penalties_scheme, ~] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
    [D_coarse, penalties_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);
    rvOperators = {}
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