The widespread use of numerical modeling codes in structural geology, particularly in forward modeling of crustal deformation, has brought valuable insight into many processes such as fault initiation and propagation, fold growth, and orogenic wedge development. Whether commercial, open source, or the result of an individual researcher’s endeavor, such numerical codes increasingly rely on advanced compiler optimizations and/or parallelization techniques. This is done to maximize performance on today’s multicore processors and thus achieve tractable model runtimes. However, such speed and productivity can come at a price. Many optimization and parallelization techniques used do indeed produce remarkable runtime speedups, but at the expense of determinism and reproducibility. Most researchers are only tangentially aware of these issues. I have aimed to explain why numerical model runs with identical sets of initial and boundary conditions can potentially produce different results and thus, sensu stricto, be nondeterministic. I have developed some solutions and guidelines for maintaining determinism and evaluated these issues using a discrete element code of upper crustal deformation. Implications for the interpretation and use of such model results in structural geology were also evaluated.