Thesis (Ph.D., Computer Science) -- University of Idaho, 2017 | Phase field modeling (PFM) is a well-known technique for simulating microstructural evolution. To model grain growth using PFM, typically each grain feature is assigned a unique non-conserved spatial variable known as an order parameter. Each order parameter field is then evolved in time. Traditional approaches for modeling these individual grains uses a one-to-one mapping of grains to order parameters since the interactions among grains is not known a priori. This presents a challenge when modeling large numbers of grains due to the computational expense of using many order parameters. This problem is exacerbated when using common numerical solution schemes including the fully-implicit finite element method (FEM), as the global matrix size is proportional to the number of order parameters squared. While previous work has developed methods to reduce the number of required variables and thus the computational complexity, none of the existing approaches can be applied to an implicit FEM implementation of PFM. Additionally, polycrystal modeling with grain growth and other coupled physics requires careful tracking of each grain's position and orientation, which is lost when using a reduced number of variables. Here, we present a modular, scalable distributed feature tracking and remapping algorithm suitable for solving these deficiencies. The algorithm presented in this dissertation maintains a unique ID for each grain even after variable remapping without restricting the underlying modeling method. This approach enables fully-coupled multiphysics using a fully generalized finite element method. Implementation details and comparative results of using this approach are presented.