This article presents a parallel algorithm to model the nonlinear dynamic interactions between ultrasonic guided waves and fatigue cracks. The Local Interaction Simulation Approach (LISA) is further developed to capture the contact-impact clapping phenomena during the wave crack interactions based on the penalty method. Initial opening and closure distributions are considered to approximate the 3-D rough crack microscopic features. A Coulomb friction model is integrated to capture the stick-slip contact motions between the crack surfaces. The LISA procedure is parallelized via the Compute Unified Device Architecture (CUDA), which enables parallel computing on powerful graphic cards. The explicit contact formulation, the parallel algorithm, as well as the GPU-based implementation facilitate LISA's high computational efficiency over the conventional finite element method (FEM). This article starts with the theoretical formulation and numerical implementation of the proposed algorithm, followed by the solution behavior study and numerical verification against a commercial finite element code. Numerical case studies are conducted on Lamb wave interactions with fatigue cracks. Several nonlinear ultrasonic phenomena are addressed. The classical nonlinear higher harmonic and DC response are successfully captured. The nonlinear mode conversion at a through-thickness and a half-thickness fatigue crack is investigated. Threshold behaviors, induced by initial openings and closures of rough crack surfaces, are depicted by the proposed contact LISA model.