def METHOD_NAME(self):
    # Randomly perturbed grid, with random linear pressure field
    g, perm = setup_cart_2d(np.array([10, 10]))
    dx = 1
    pert = 0.4
    g.nodes = g.nodes + dx * pert * (
        0.5 - np.random.rand(g.nodes.shape[0], g.num_nodes)
    )
    # Cancel perturbations in z-coordinate.
    g.nodes[2, :] = 0
    g.compute_geometry()
    bound_faces = np.argwhere(np.abs(g.cell_faces).sum(axis=1).A.ravel("F") == 1)
    bound = pp.BoundaryCondition(
        g, bound_faces.ravel("F"), ["dir"] * bound_faces.size
    )
    # Python inverter is most efficient for small problems
    flux, bound_flux, *_ = pp.Mpfa("flow")._flux_discretization(
        g, perm, bound, inverter="python"
    )
    div = g.cell_faces.T
    a = div * flux
    pr_bound, pr_cell, gx, gy = setup_random_pressure_field(g)
    rhs = div * bound_flux * pr_bound
    pr = np.linalg.solve(a.todense(), -rhs)
    p_diff = pr - pr_cell
    self.assertTrue(np.max(np.abs(p_diff)) < 1e-8)