function diffOpsStruct = constructDiffOpsMultiStage(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_unstable = 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_unstable', D_unstable,...
                     'D_flux', D_flux,...
                     'D_t', D_t,...
                     'penalties_scheme', penalties_scheme,...
                     'penalties_res', penalties_res);
end