  At mesoscopic scales electrolyte solutions are modeled by the fluctuating
generalized Poisson-Nernst-Planck (PNP) equations [J.-P. P\'eraud et al., Phys.
Rev. F, 1(7):074103, 2016]. However, at length and time scales larger than the
Debye scales, electrolytes are effectively electroneutral, and the
charged-fluid PNP equations become too stiff to solve numerically. Here we
formulate the isothermal incompressible equations of fluctuating hydrodynamics
for reactive multispecies mixtures involving charged species in the
electroneutral limit, and design a numerical algorithm to solve these
equations. Our model does not assume a dilute electrolyte solution but rather
treats all species on an equal footing, accounting for cross-diffusion and
non-ideality using Maxwell-Stefan theory. By enforcing local electroneutrality
as a constraint, we obtain an elliptic equation for the electric potential that
replaces the Poisson equation in the fluctuating PNP equations. We develop a
second-order midpoint predictor-corrector algorithm to solve either the
charged-fluid or electroneutral equations with only a change of the elliptic
solver. We use the electroneutral algorithm to study a gravitational fingering
instability, triggered by thermal fluctuations, at an interface where an acid
and base react to neutralize each other. Our results demonstrate that, because
the four ions diffuse with very different coefficients, one must treat each ion
as an individual species, and cannot treat the acid, base, and salt as neutral
species. This emphasizes the differences between electrodiffusion and classical
Fickian diffusion, even at electroneutral scales.
