This paper gives an overview of some recent works by the author and collaborators on numerical methods for partial differential equations. One focus is on the development and analysis of algebraic solvers (especially multigrid methods) for discretized partial differential equations(PDEs) by closely taking into account qualitative properties of PDEs and their finite element discretizations together with grid adaptation. We will give several examples to demonstrate that optimal PDE solvers can be designed to be user-friendly and practical, robust with respect to various physical parameters and applicable to a variety of linear and nonlinear problems. More specifically, we will present a sharp convergence theory for the method of subspace corrections (such as multigrid and domaindecomposition) for SPD systems. In particular, we demonstrate that the method of subspace correction will be robust for nearly singular systems if the space decomposition satisfies the near-null space recovery condition (namely, the near-null space can be fully recovered by the parts of near-null space contained in all subspaces). We present a solution of an open problem on robust multigrid and domain decomposition methods for problems with strongly discontinuous jumps. We discuss a new class of optimal and practical (AMG) preconditioners for H(curl) and H(div) systems. We present a methodology that reduces the solution of an indefinite system to one or very few (nearly singular) SPD systems (which can be solved efficiently by the method of subspace correction). Finally, we use a class of non-Newtonian models to illustrate how to design special discretization schemes that are solver-friendly and also satisfy other desirable properties.