# Should I expect computational gains using a second-order splitting method here?

I am trying to solve a three-dimensional baroclinic transport problem. The hydrodynamic (three-dimensional shallow water) equations are: \begin{align} \nabla \cdot \vec{v} = 0, \tag{1} \\ \frac{\partial u}{\partial t} + \nabla \cdot (u \vec{v} ) + \frac{1}{\rho_0} \frac{\partial P}{\partial x} - f_{\text{coriolis}}v - \frac{1}{\rho_0}(\nabla \cdot T_{\text{mx}}) - S_x = 0, \tag{2} \\ \frac{\partial v}{\partial t} + \nabla \cdot (v \vec{v} ) + \frac{1}{\rho_0} \frac{\partial P}{\partial y} + f_{\text{coriolis}}u - \frac{1}{\rho_0}(\nabla \cdot T_{\text{my}}) - S_y = 0, \tag{3} \\ \frac{\partial H}{\partial t} + \frac{\partial (\bar{u}H)}{\partial x} + \frac{\partial (\bar{v}H)}{\partial y} = 0 \tag{4}. \end{align} The constituent (salinity or temperature) transport equation is: \begin{align*} \frac{\partial c}{\partial t} = \nabla \cdot (\vec{v} c) + \nabla \cdot (\bar{\bar{D}} \nabla c) = 0 \tag{5}. \end{align*} The problem is to solve for water velocity $\vec{v} = [u,v,w]^T$, river depth $H$, and constituent concentrate $c$ over a fluid domain $\Omega \subset \mathbb{R}^3$ and timescale $[0, t_{\text{end}}]$ of interest. The first equation is usually called the incompressible continuity equation. The second and third equations are called the horizontal momentum equations. The fourth equation is called the depth-averaged continuity equation. The fifth equation is the transport equation.

The other important relationship to note is that the fluid pressure $P$ (which appears in the horizontal momentum equations) is related to the concentration $c$. Specifically, $P$ depends on the fluid density $\rho$, and $\rho$ is related to $c$ through an algebraic equation of state. So there is two-way coupling between the hydrodynamic equations and the transport equation.

The solution procedure used to run the simulation one time-step is as follows. Note that spatial discretization is done using a first-order SUPG FEM, and temporal discretization is done using a second-order (implicit) BDF method.

1. Using $u,v,w,H,P$ from the old time-step, we solve discrete weak version of the horizontal momentum equations (2)-(3) and the depth average continuity (4) monolithically to get $u, v, H$ at the new time-step. Since the horizontal momentum equations (2)-(3) are nonlinear, Newton's method with line-search is used.
2. Using the new $u, v, H$ from step 1, a discrete weak version of incompressible continuity equation (1) is solved to get $w$ at the new time-step. This is a linear problem, so this is computationally easy compared to step 1.
3. Using the new $u, v, w, H$ from steps 1 and 2, a discrete weak version of the transport equation (5) is solved to get $c$ at the new time-step. This is also linear, so this is computationally easy compared to step 1. Also, this new $c$ value is used in the calculation of a new $P$ value for use in step 1.

After these three steps, I have all the variables we are looking for at the next discrete time value for all spatial coordinates in our grid. So I repeat these three steps to march in time.

This solution procedure might be seen as a first-order operator splitting schemes, since the hydrodynamics processes (equations (1)-(4) and steps 1-2) are solved separately from the transport process (equation (5) and step 3). Should one expect to gain any sort of advantage (e.g., in computational efficiency) by implementing second-order splitting schemes (e.g., from this webpage: http://www.asc.tuwien.ac.at/~winfried/splitting/)?

I have ran some computational experiments testing a few second-order schemes where the hydrodynamic equations (1)-(4) are split from equation (5), i.e., the coefficients $a_i$ correspond fractional step coefficients for hydrodynamic processes while the coefficients $b_i$ correspond to fractional step coefficients for transport processes; however, I do not seem to see any computational gains at the moment.

The experiments are run on a lock-exchange problem that is similar to the one found in [Shin et al. 2014]. Initially, the left-half of the flume (with length 2 meters, width 0.2 meters, and depth 0.2 meters) is filled with saltwater, the right-half is filled with freshwater, so there is a sharp interface in the middle at $t = 0$. The barrier is instantly removed, so the dense saltwater flows beneath the less-dense freshwater. The fixed uniform mesh that I use to run my experiments is relatively coarse (grid spacing is 0.2 meters), but I don't think this matter since I am only interested in temporal convergence. I roughly measure computationally efficiency by plotting the discrete $L_2$ error versus the wall-clock time of the simulation. Below is an example of one such plot (split1 through split5 correspond to different implementations of second-order splitting schemes, while split0 corresponds to first-order splitting). Note that the "reference solution" used in this discrete $L_2$ error measurement is actually an approximate solution to the semi-discrete system, which is computed using a very small time-step size of 0.01 seconds. Also, I am specifically looking at the error in salt concentration $c$. Looking at points which are vertically aligned, you can see that the black points (first-order splitting) achieves an $L_2$ error which is roughly the same as the $L_2$ errors computed using second-order splitting (in fact, they are a little better).

References:

Shin, J. O., S. B. Dalzieland, and P. F. Linden. 2004. Gravity currents produced by lock exchange. Journal of Fluid Mechanics 521:1–34.