Source(src,double(c2),-rate,'sat',[1,1]);
    
    %% Set up solution structures with only one phase
    refSol = initResSol(g_fine, bhp, [0, 1]);
    msSol = refSol;

    gravity off
    verbose = false;
    
    %% Set up pressure and transport solvers
    
    % Reference TPFA
    r_psolve = @(state) incompTPFA_DFM(state, g_fine, S, fluid, 'Wells', W,'src',src,'c2cTrans',S2,'MatrixOutput',true);
    
    
    % MsFV using a few GMRES iterations to improve flux error
    psolve   = @(state) solveMSFV_TPFA