An adjoint Petrov-Galerkin method is proposed which solves multidimensional advection-dispersion, diffusion, and Poisson equations by means of weight (test) functions that form numerical solutions of adjoint state equations on a sequence of nested grids. The cells of these grids may have straight or curved boundaries. The weights are numerical relatives of 'optimal test functions' derived analytically by Celia et al. but possess important additional attributes. Solutions corresponding to a given member of the sequence are ideally suited for parallel computation. All solutions with the possible exception of those corresponding to the finest local grids are free of advective terms and contain diffusive terms only along intercell boundaries, a property which enhances accuracy and makes the method especially well suited for the treatment of nonuniform velocity fields in heterogeneous media. While the dependent variable is evaluated only at the nodes of the coarset grid, its spatial variation is described with a resolution corresponding to the finest grids. This is achieved in part by associating with each member of the sequence a novel family of 'optimal interpolation (trial) functions' which depend solely on the corresponding weight functions and reflect the general asymmetry of the latter. Both the weight and interpolation functions can vary with time to maintain optimality when the velocity and parameter fields are not constant.