2.4.4. Laplacian Transformation

*f*

GCCOM solves the full nonhydrostatic 3D Navier–Stokes equations as follows:

$$\frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} = -\frac{1}{\rho\_0} \nabla p - \frac{\mathcal{g}\rho}{\rho\_0} \vec{k} - \nabla \cdot \vec{\overline{\tau}},\tag{6}$$

$$\frac{\partial T}{\partial t} + \vec{u} \cdot \nabla T = \nabla \cdot (k\_T \nabla T),\tag{7}$$

$$\frac{\partial \mathcal{S}}{\partial t} + \vec{u} \cdot \nabla \mathcal{S} = \nabla \cdot (k\_{\mathcal{S}} \nabla \mathcal{S}),\tag{8}$$

$$
\nabla \cdot \vec{u} = 0,
$$

$$
\rho = \rho(T, S, p). \tag{10}
$$

Here, *u* = (*<sup>u</sup>*, *v*, *w*) are velocities, *gρρ*0*k* is gravity acceleration, -*τ* our stress tensor solved by Large Eddie Simulation, *ρ* is a equation of state and *kT*,*<sup>S</sup>* diffusivity constants for temperature and salinity.

The Boussinesq approximation [48] is applied to Equation (6) by taking the divergence and substituting ∇*u* as in Equation (9). This step cancels every term other than the pressure *p* from Equation (6), leaving a homogeneous Laplacian problem to solve as in Equation (11), which by itself is computationally expensive to solve, as we will discuss in the rest of this section:

$$
\nabla^2 p = \frac{\partial^2 p}{\partial x} + \frac{\partial^2 p}{\partial y} + \frac{\partial^2 p}{\partial z} = 0,
\tag{11}
$$

$$\begin{split} \nabla^2 p &= L(p) - L(\mathbf{x}) \left[ \xi\_x \frac{\partial p}{\partial \xi} + \eta\_x \frac{\partial p}{\partial \eta} + \zeta\_x \frac{\partial p}{\partial \zeta} \right] \\ &\quad - L(\mathbf{y}) \left[ \xi\_y \frac{\partial p}{\partial \xi} + \eta\_y \frac{\partial p}{\partial \eta} + \zeta\_y \frac{\partial p}{\partial \zeta} \right] \\ &\quad - L(\mathbf{z}) \left[ \zeta\_z \frac{\partial p}{\partial \zeta} + \eta\_z \frac{\partial p}{\partial \eta} + \zeta\_z \frac{\partial p}{\partial \zeta} \right]. \end{split} \tag{12}$$

The GCCOM Laplacian in curvilinear coordinates is formulated in [29]; note that we are following notation from [31], the operator *L*() is defined in Equation (13), where *a*, *b*, *c*, *d*,*e*, *q* are coefficients related to the curvilinear transformation [12] that are defined in terms of the Jacobian *J* and generalized in two sets of coordinates by their even commutation: {*<sup>x</sup>*, *y*, *z*} as {1, 2, <sup>3</sup>}, and {*ξ*, *η*, *ζ*} as {*<sup>A</sup>*, *B*, *<sup>C</sup>*}, a general rule for the derivatives can be defined as in Equation (14), where we have adopted a compact differential representation, i.e., *A*1 = *∂ξ∂x* ,

$$L(\cdot) = a\frac{\partial^2(\cdot)}{\partial \xi^2} + b\frac{\partial^2(\cdot)}{\partial \eta^2} + c\frac{\partial^2(\cdot)}{\partial \zeta^2} + 2\left[d\frac{\partial^2(\cdot)}{\partial \xi \partial \eta} + e\frac{\partial^2(\cdot)}{\partial \zeta \partial \eta} + q\frac{\partial^2(\cdot)}{\partial \xi^2 \partial \zeta}\right],\tag{13}$$

$$A\_1 = \frac{\partial \mathfrak{J}}{\partial \mathbf{x}} = \mathfrak{J}(\mathfrak{L}\_B \mathfrak{J}\_\mathbb{C} - \mathfrak{L}\_\mathbb{C} \mathfrak{J}\_B). \tag{14}$$

By applying these elements to the transformed Laplacian operator and discretizing using 2nd-order finite difference, Equation (15) becomes the discretized transformed Laplacian operator in general curvilinear coordinates [33],

$$\begin{split} \nabla^2 p &= -\frac{1}{2(\Delta\_{\xi}^2 \Delta\_{\eta}^2 \Delta\_{\xi}^2)} \Bigg( \mathbf{A} \mathbf{z}(i,j,k) \boldsymbol{p}(i,j,k) \\ &+ [\beta\_1(i,j,k) + \beta\_2(i,j,k)] \boldsymbol{p}(i+1,j,k) \\ &+ [\beta\_1(i,j,k) - \beta\_2(i,j,k)] \boldsymbol{p}(i-1,j,k) \\ &+ [\lambda\_1(i,j,k) + \lambda\_2(i,j,k)] \boldsymbol{p}(i,j+1,k) \\ &+ [\lambda\_1(i,j,k) - \lambda\_2(i,j,k)] \boldsymbol{p}(i,j-1,k) \\ &+ [\tau\_1(i,j,k) + \tau\_2(i,j,k)] \boldsymbol{p}(i,j,k+1) \\ &+ [\tau\_1(i,j,k) - \tau\_2(i,j,k)] \boldsymbol{p}(i,j,k-1) \\ &- \tau\_{\text{xy}}(i,j,k) \big( p(i+1,j+1,k) + p(i-1,j-1,k) \big) \\ &+ \tau\_{\text{xy}}(i,j,k) \big( p(i+1,j-1,k) + p(i-1,j+1,k) \big) \\ &- \tau\_{\text{yz}}(i,j,k) \big( p(i,j-1,k-1) + p(i,j+1,k+1) \big) \\ &+ \tau\_{\text{yz}}(i,j,k) \big( p(i,j-1,k+1) + p(i,j+1,k-1) \big) \\ &- \tau\_{\text{xz}}(i,j,k) \big( p(i-1,j,k-1) + p(i+1,j,k+1) \big) \\ &+ \tau\_{\text{xz}}(i,j,k) \big( p(i-1,j,k+1) + p(i+1,j,k-1) \big) \Biggr), \\ \end{split}$$

where *α*, *β*1, *β*2, *λ*1, *λ*2, *τ*1, *τ*2, *<sup>τ</sup>xy*, *<sup>τ</sup>yz*, *τxz* are transformation coefficients found after algebraic manipulation.

Now, the pressure can be solved entirely as a system of linear equations in the form *Ax* = *b*, where *A* is the system matrix constructed by the position coefficients, *b* is the known pressure on each point (Right Hand Side, or *RHS*), and *x* is the solution vector for pressure.

A peculiarity of this system is the inherent lexicographical ordering: the points of the 19-point stencil follow a specific positioning pattern in the program memory, as seen in Figure 2: on each *xy*-plane, points are sorted in the *y*-direction before the *x*-direction, then the points are sorted by planes on the bottom *z*-axis first; this is the same as a front-to-back, bottom-up, ordering in the 3D stencil. This bookkeeping is crucial to obtain the right dynamics out of the Laplacian.

**Figure 2.** Lexicograpical ordering of the 19-point stencil. Here, element 10 is the current point of the stencil [38].

It is important to note that the coefficient matrix, *A* is large [*gx*, *gy*, *gz*] × [*gx*, *gy*, *gz*] , sparse, and in general not singular, and non-symmetric [33]. This automatically eliminates methods that directly invert the matrix to solve *x* = *A*−1*b* for large problem sizes. In GCCOM, this equation is solved with the Generalized Conjugate Residual (GCR) method preconditioned with block Jacobi. At this point, we are able to take advantage of the PETSc linear solver system which includes a long list of solvers and preconditioners that can be accessed via command line arguments: -ksp\_type [solver] -pc\_type [precond].

#### 2.4.5. External Boundary Data

External boundary conditions are generated outside of GCCOM and stored separately in ASCII files. They contain velocity data at the boundaries for planes in the *x*-direction (west and east) and *z*-direction (north and south). Processors located on these boundaries read boundary condition data into local memory. As we work with higher resolution meshes, this external file reading becomes an obstacle to overcome and adds to the I/O overhead. Part of the required future work includes updating these routines to read parallel NetCDF files.

#### *2.5. Test Case Experiments*

In this section, we describe the experiments used to validate the parallel framework of GCCOM. The Seamount experiment was chosen for being the most comprehensive for the model, using most—if not all—of the model capabilities at once.
