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$.

enter image description here 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).


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


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.