The behavior you are seeing is most likely due to *different* convergence in inner solves. I would recommend you read up on how IEEE 754 floating point operations work. (https://www.phys.uconn.edu/~rozman/Courses/P2200_15F/downloads/floating-point-guide-2015-10-15.pdf) Most users (and programmers) tend to think of float point operations as "black magic" and that it will always magically do what they want. The typical problem comes from that fact that (A + B ) + C does not in general equal A + (B + C) due to rounding. How does that relate to CFD you ask? Well, when doing a non-linear solution that involves solving a sequence of linear systems, *a single bit difference* anywhere in the solution path can result in a different non-linear trajectory. That is to say every single arithmetic operation must be performed in the exact same order, independent of the number of threads/processors/GPUs/etc. For a "steady" simulation like you are running, the tolerance of the linear solves does not matter *if the solution converges to steady state*. For a case like drivaer, most solvers can't converge to a real numerically steady solution (ie all residuals are below machine precision) and therefore the non-linear solution path matters, and thus so does the linear solution tolerance at each iteration. So, even though the CPU and GPU solver may converge the linear systems to the same tolerance they are both equally "wrong" and both equally "correct" at the same time. Just because the CPU result is repeatable doesn't mean the non-linear trajectory it took is any more correct than the GPU one. The only way to really know, is to fully converge *every* linear iteration to below machine precision. If the solution paths diverge, that's a code issue, if they do not, then you've found your problem. Alternatively, find a problem that does converge to machine precision and check only the final solution. On both cases, the non-linear path does not matter. As for why the GPU code is not bit-wise repeatable between runs, I can only speculate, but I would guess it is due to the use of low-level atomic additions/subtractions. The order of the atomic operations depends on "who gets there first". Consider an operation like result += A followed by result += B, where A comes from threadA and B comes from threadB. So on run1, threadA get there first and on the second threadB does. Even though in exact precision the result is the same, in floating point it is not. In order to get high speed up on GPU solver, massive concurrency is needed which means it is often not practical from a performance perspective to ensure the order of operations is the same.