Numerically Solving Fluid Dynamics I'm trying to model the behavior of a scalar field $\rho(\vec x)$.  Associated with $\rho$ there is an energy density $\mu$ like so:
$$
\mu = \tau[\rho,Q] + u[\rho,Q]
$$
$$
u = u_{ext}[\rho,Q] + u_{ee}[\rho]
$$
The symbols $\tau$, $u$, $u_{ext}$ and $u_{ee}$ represent scalar fields that will be computed from $\rho$ and the knowledge of a set of "point charges" $Q$ that don't change over time but which influence the dynamics of $\rho$. The terms between square brackets represent the dependency of each field in the equations with either $\rho$ or $Q$.
I want to study two types of equations of motion, the first one is a continuum version of Newton's second law of motion:
$$
\partial_t \vec v = -\nabla \mu
$$
The second type is
$$
\partial_t \rho = -\nabla \mu
$$
I understand that these two equations describe different types of motion and one of the points of my study is to determine which one is more appropriate for my purposes, so ideally I would like to program both.
In order to do that I have to know how to compute $\rho_{t+\delta t}$ as a function of $\rho_{t}$ from either of the equations above. I read somewhere that in order to determine the dynamics of a fluid I have to solve a set of differential equations called the Navier-Stokes Equations.
I understand that to describe the dynamics of a set of point particles following classical mechanics, for example, formally I have to solve Newton's equations of motion. But when you actually sit down to write a software that does that, there are actually a number of algorithms you can use, such as Verlet or Leap-frog, that generate trajectories which are solutions to Newton's equations. I was wondering whether there was any such algorithm I could similarly implement to solve my hydrodynamical problem.
So my question is:
What time-stepping algorithms are available to model the dynamics of a 3D field like this one?
 A: While hydrodynamic equations are represented as a coupled set of PDEs, the time stepping is generally treated as an ODE instead,
$$\frac{\mathrm d\mathbf u}{\mathrm d t}=F(\mathbf u(t))$$
where $\mathbf u=(\rho,\,\pi_x,\,\pi_y,\,\pi_z,\,E)^T$ is the state vector ($\rho$ the density, $\pi_i=\rho v_i$ the momentum density in the direction $i$ and $E$ the total energy) and $F(\cdot)$ the flux. Hence, any of the methods used in solving an ODE are applicable to hydrodynamics.
Some of the more common increments I have seen in academic codes are,

*

*Forward Euler method update
$$\rho(t+\Delta t)=\rho(t)-\nabla\mu(t)\cdot\Delta t$$
which is ubiquitous due to its simplicity.

*

*Backward Euler method,
$$\rho(t+\Delta t)=\rho(t)-\nabla\mu(t+\Delta t)\cdot\Delta t,$$
typically involves matrix operations/reductions and are thus more costly than the forward update, but is still an option.



*Predictor-corrector method, which solves the state at $t+\Delta t$ using a predicted value & an interpolated value,
\begin{align}
  \hat{\rho}(t+\Delta t)&=\rho(t)-\nabla\mu(t,\,\rho(t))\cdot\Delta t \\
  \rho(t+\Delta t)&=\rho(t)-\frac{1}{2}\left(\nabla\mu(t+\Delta t,\,\hat{\rho})+\nabla\mu(t,\,\rho)\right)\cdot\Delta t
\end{align}

*Linear multistep methods, most notably  Runge-Kutta methods, are pretty common as well, though in practice I have not seen anything beyond RK2 used in academic hydrodynamic codes. Note that the more steps one includes in these multistep methods, the longer each step of the simulation will take.

Note that, in hydrodynamics, your time step $\Delta t$ is constrained by the Courant-Friedrichs-Lewy condition, which essentially argues that the rate of information transfer (e.g., speed of sound) limits how far material can travel from one cell to another in one time step, lest unphysical negative densities/pressures or superluminal velocities occur. This conditions enforces that,
$$\Delta t\leq \mathcal{C}\frac{\Delta x}{u_\text{max}}$$
where $u_\text{max}$ is the maximum characteristic speed over the whole domain at the current point and $\mathcal{C}\leq1$ is a constant for the whole simulation. One can derive such time step constraints by considering the von Neumann stability analysis.
