We present Rivus, a differentiable two-dimensional incompressible Navier-Stokes solver built on JAX and XLA that enables gradient-based inverse design of fluid flow through automatic differentiation (AD) of the full time integration. The solver combines a pseudo-spectral pressure projection with semi-Lagrangian BFECC advection in a Strang splitting framework, achieving unconditional stability and second-order temporal accuracy—properties essential for stable differentiation through long trajectories. We introduce a stream function parameterization that enforces the divergence-free constraint by construction, eliminating the need for projection in the optimization path, and employ Sobolev gradient preconditioning with spectral annealing for stable convergence. A multiscale continuation method provides an asymptotically coarsened optimization landscape with wider basins of attraction. We demonstrate the method on two problems: (i) optimal velocity field design for density transport, and (ii) vortex-optimized mixing, where we establish a connection to Prigogine’s theory of dissipative structures. In the mixing problem, the optimizer spontaneously discovers organized vortex topologies—from simple two-roll states to multi-scale filamentary cascades—as the optimal route to entropy production, with sharp bifurcations in the optimal control space as viscosity and simulation horizon are varied. Additional considerations necessary for the stable propagation of gradients through large computational graphs representative of 10,000+ timesteps will be elucidated. The solver is parallelized across multiple AMD Instinct GPUs via mesh sharding and achieves high performance and efficient scaling on the MI300A and MI355X architectures. Finally, we discuss extensions to learned turbulence closures via differentiable super-resolution and future studies.