  Partial differential equations (PDEs) have become an essential tool for
modeling complex physical systems. Such equations are typically solved
numerically via mesh-based methods, such as finite element methods, with
solutions over the spatial domain. However, obtaining these solutions are often
prohibitively costly, limiting the feasibility of exploring parameters in PDEs.
In this paper, we propose an efficient emulator that simultaneously predicts
the solutions over the spatial domain, with theoretical justification of its
uncertainty quantification. The novelty of the proposed method lies in the
incorporation of the mesh node coordinates into the statistical model. In
particular, the proposed method segments the mesh nodes into multiple clusters
via a Dirichlet process prior and fits Gaussian process models with the same
hyperparameters in each of them. Most importantly, by revealing the underlying
clustering structures, the proposed method can provide valuable insights into
qualitative features of the resulting dynamics that can be used to guide
further investigations. Real examples are demonstrated to show that our
proposed method has smaller prediction errors than its main competitors, with
competitive computation time, and identifies interesting clusters of mesh nodes
that possess physical significance, such as satisfying boundary conditions. An
R package for the proposed methodology is provided in an open repository.
