  In this paper we propose an extension of the generalized Lagrangian
multiplier method (GLM) of Munz et al. (JCP 2000, JCP 2002), which was
originally conceived for the numerical solution of the Maxwell and MHD
equations with divergence-type involutions, to the case of hyperbolic PDE
systems with curl-type involutions. The key idea here is to solve an augmented
PDE system, in which curl errors propagate away via a Maxwell-type evolution
system. The new approach is first presented on a simple model problem, in order
to explain the basic ideas. Subsequently, we apply it to a strongly hyperbolic
first order reduction of the CCZ4 formulation (FO-CCZ4) of the Einstein field
equations of general relativity, which is endowed with 11 curl constraints.
Several numerical examples, including the long-time evolution of a stable
neutron star in anti-Cowling approximation, are presented in order to show the
obtained improvements with respect to the standard formulation without special
treatment of the curl involution constraints. The main advantages of the
proposed GLM approach are its complete independence of the underlying numerical
scheme and grid topology and its easy implementation into existing computer
codes. However, this flexibility comes at the price of needing to add for each
curl involution one additional 3 vector plus another scalar in the augmented
system for homogeneous curl constraints, and even two additional scalars for
non-homogeneous curl involutions. For the FO-CCZ4 system with 11 homogeneous
curl involutions, this means that additional 44 evolution quantities need to be
added.
