Tide-driven flow at Bight of Abaco (Bahamas)¶
In this example, we examine the simulation of a tide-driven using a quadrature-free discontinuous Galerkin (DG) discretization of the 2D shallow water equations (SWE).
Problem description¶
We consider the 2D shallow water equations in conservative form defined on some 2D domain \(\Omega\).
Denoting \(\boldsymbol{c} := (\xi, U,V)^T\), the system is given in the following compact form:
where
This example uses the following types of boundary conditions for the SWE:
Land: No normal flow \(\boldsymbol{q} \cdot \boldsymbol{n} = 0\).
Open-sea: Prescribed free surface elevation \(\xi(t, \boldsymbol{x}) = \hat{\xi}(t, \boldsymbol{x})\).
The following boundary conditions are furthermore available for different scenarios:
Dirichlet: All unknowns are specified \(\boldsymbol{q}(t, \boldsymbol{x}) = \hat{\boldsymbol{q}}(t, \boldsymbol{x}), \;\; \xi(t, \boldsymbol{x}) = \hat{\xi}(t, \boldsymbol{x})\).
River: For supercritical flow examples, we set the following river (inflow) boundary conditions: \(q_{\boldsymbol{n}}(t, \boldsymbol{x}) = \hat{q}_{\boldsymbol{n}}(t, \boldsymbol{x}), \quad q_{\boldsymbol{\tau}}(t, \boldsymbol{x}) = \hat{q}_{\boldsymbol{\tau}}(t, \boldsymbol{x}), \quad \xi(t, \boldsymbol{x}) = \hat{\xi}(t, \boldsymbol{x})\) with the normal and tangential integrated velocities \(\hat{q}_{\boldsymbol{n}}(t, \boldsymbol{x})\) and \(\hat{q}_{\boldsymbol{\tau}}(t, \boldsymbol{x})\).
Radiation: At the outflow boundary of supercritical flow examples, we specify radiation boundary conditions where no unknowns are prescribed.
The initial conditions for the elevation and integrated velocity can be provided as follows:
\(\xi(0, \boldsymbol{x}) = \xi_0(\boldsymbol{x}),\quad {\boldsymbol{q}}(0, \boldsymbol{x}) = {\boldsymbol{q}}_0(\boldsymbol{x})\)
This example uses the lake-at-rest initial condition.
Quadrature-free formulation¶
To avoid fracture-type nonlinearities to permit a quadrature-free DG formulation, we modify the system as follows (cf. [FKAZGK2020]). Introducing the depth averaged velocity \(\boldsymbol{u}=(u,v)^T\) related to \(\boldsymbol{q}\) by \(\boldsymbol{q} = \boldsymbol{u}H\) the new system reads
where \(H^2-h_b^2=\xi^2+2h_b\xi+h_b^2-h_b^2=\xi(H+h_b)\) and thus
DG discretization¶
As a next step, we discretize the continuous problem. Correspondingly, the physical domain \(\Omega\) is triangulated into triangles \(\Omega_e, e \in \{0,...,E\}\). The local variational formulation of the system on an element \(\Omega_e\) is obtained by multiplying with sufficiently smooth test functions \(\boldsymbol{\phi}\) and \(\boldsymbol{\psi}\) followed by the integration by parts where we use the notation \((\cdot,\cdot)_{\Omega_e}\) and \(\langle\cdot,\cdot\rangle_{\partial \Omega_e}\) for the \(L^2\)-scalar products on elements and edges, respectively, and denote by \(\boldsymbol{n}_e\) an exterior unit normal to \(\partial \Omega_e\)
Next for obtaining the semi-discrete formulation, let \(\mathbb{P}^p(\Omega_e)\) be the polynomial space of order \(p\) on \(\Omega_e\). We replacing \(\boldsymbol{c}\) and \(\boldsymbol{u}\) by their finite-dimensional counterparts \(\boldsymbol{c}_\Delta, \boldsymbol{u}_\Delta\) and utilize test functions \(\boldsymbol{\phi}_\Delta \in \mathbb{P}^p(\Omega_e)^3, \, \boldsymbol{\psi}_\Delta \in \mathbb{P}^p(\Omega_e)^2\):
The edge flux \(\boldsymbol{{A}}(\boldsymbol{c}_\Delta, \boldsymbol{u}_\Delta) \cdot\boldsymbol{n}_e\) is approximated on \(\partial\Omega_e\) by a numerical flux \(\boldsymbol{\widehat{A}}(\boldsymbol{c}_\Delta, \boldsymbol{u}_\Delta,\boldsymbol{\tilde{c}}_\Delta, \boldsymbol{\tilde{u}}_\Delta, \boldsymbol{n}_e)\) that depends on discontinuous values of the solution on element \(\Omega_e\) (without tilde) and on its edge neighbours (with tilde). On exterior domain boundaries, the specified elevation and velocity boundary values are utilized in the flux computation. In this example, we use the Lax–Friedrichs flux, defined as
where \(\lambda\) denotes the largest (in absolute value) eigenvalue of \(\frac{\partial}{\partial \boldsymbol{c}} \boldsymbol{\widetilde{A}}(\boldsymbol{c}) \cdot \boldsymbol{n}_e\). Denoting by \(\boldsymbol{x}_f\) the midpoint of edge \(f \subset \partial\Omega_e\), our quadrature-free scheme uses the following approximation of \(\lambda|_f\) (cf. [FKAZGK2020]):
We now take a basis \(\varphi_{ei}(\boldsymbol{x}), \,i=1, \ldots, K(p)\), of \(\mathbb{P}^p(\Omega_e)\) and can represent \(\boldsymbol{c}_\Delta\) and \(\boldsymbol{u}_\Delta\) as
with \(\boldsymbol{e}_j\) denoting the \(j\)-th unit vector in \(\mathbb{R}^3\) or \(\mathbb{R}^2\). A basis of \(\mathbb{P}^p(\Omega_e)\) is defined via a mapping from the corresponding reference basis on \(\hat{\Omega}\)
where \(F_e: \hat{\Omega}\rightarrow\Omega_e,\; \hat{\boldsymbol{x}}\rightarrow \boldsymbol{x} := \boldsymbol{B}_e\hat{\boldsymbol{x}}+\boldsymbol{a}_{e1}\) is the affine-linear transformation from the reference triangle onto triangle \(\Omega_e\) as shown in the figure below with
and \(\boldsymbol{a}_{e1}, \boldsymbol{a}_{e2}, \boldsymbol{a}_{e3}\) denoting the vertex coordinates of \(\Omega_e\).
The number of basis functions \(K(p)\) is dependent on the respective polynomial space and has the following values in \(\mathbb{R}^2\): \(K(0)\!=\!1, \, K(1)\!=\!3\) and \(K(2)\!=\!6\). The reference basis functions employed in our implementation are orthonormal with respect to the \(L^2\)-scalar product on \(\hat{\Omega}\) and given by
Up to this point, we defined the discrete formulation of the problem we want to solve. The next step is to look at the algebraic representation of the element and edge integrals that arise.
For compactness, we show an algebraic representation of our discrete scheme only for element integrals. We insert the basis representations into the system and test the first equation with \(\boldsymbol{\phi}_\Delta = \varphi_{ek} \boldsymbol{e}_1\) and obtain for \(k \in {1,\dots,K(p)}\):
For \(\boldsymbol{\phi}_\Delta = \varphi_{ek} \boldsymbol{e}_2\) and \(k \in {1,\dots,K(p)}\), we get:
And for \(\boldsymbol{\phi}_\Delta = \varphi_{ek} \boldsymbol{e}_3\) we get a similar expression.
Mapping to code¶
As ExaStencils does not directly support higher-order discretizations like the one we are using, we adapted the workflow as shown below. Instead of starting with the continuous model in the multi-layered approach of ExaSlang, the continuous and discrete problems are manually derived and implemented in GHODDESS, which then directly maps to layer 4 [FKAZGK2020].
We want to examine the three steps of mapping the discretized equations to efficient simulation code.
Math to Python¶
In this step, the DG scheme is mapped to Python code with the help of the SymPy package. This is facilitated by already provided basic abstractions for, e.g., basis functions, just as classes representing triangles and data fields. Since ExaStencils is restricted to quadrilateral meshes, we divide each element into two differently oriented triangles – lower and upper – to obtain a triangular mesh for our discretization.
The algebraic representations of the discrete scheme shown above are implemented manually in our framework. Here we consider the first part of the element integral for \(\xi\) as an example:
which can be implemented in GHODDESS by
sum(cu(tri.orientation, k) * tri.det_b_inv * (
tri.b[1, 1]
* integrate_over_tri(basis.phi[k] * basis.phi_dx[p])
- tri.b[1, 0]
* integrate_over_tri(basis.phi[k] * basis.phi_dy[p]))
for k in range(d))
Information about the triangles, such as the orientation or the determinant of the mapping from the reference triangle and its reciprocal, are stored in \(\texttt{tri}\). SymPy is then used to evaluate integrals and derivatives analytically.
Python to ExaSlang¶
To allow a mapping from the symbolic representation to ExaSlang layer 4 code, SymPy expressions are enriched with a few required abstractions such as field symbols that correspond to accesses to ExaSlang fields storing quantities defined on the computational domain. GHODDESS then creates an auxiliary knowledge file holding parameters that guide the generation process, as well as ExaSlang specifications.
As an example, we once again look at the implementation of the element integral for \(\xi\) from above. Depending on the chosen order, this expression will be evaluated differently. For linear order and the second degree of freedom, i.e., \(\texttt{d=3}\) and \(\texttt{p=2}\), the underlying sum over the integrals can be evaluated analytically as
using the basis formulation from above. Together with the remaining terms, this representation can then be mapped to ExaSlang. We now consider only the first half of the element integral and ignore the second half, edge integrals, and the right-hand side. The corresponding update kernel would look like shown below.
loop over cxiNewLower1 {
cxiNewLower1 = ( -3.0 * 1.41421356237310
* bLower3 * cuLower0 * invDetBLower0 )
}
Frequently used sub-expressions, such as parts of the transformation matrix \(B\) and the reciprocal of its determinant, are stored in designated data fields. This kernel only updates triangles with lower orientation and in practice an equivalent kernel updating the remaining upper triangles is generated.
ExaSlang to code¶
The ExaStencils code generator is capable of parsing the emitted ExaSlang code, applying code transformations and certain optimizations, and, finally, outputting a C++ or CUDA code combined with MPI for distributed memory parallelism.
Following the example from the previous paragraph, we can now generate compilable C++ code similar to the one depicted below without optimizations to give developers insight.
for (int i1 = 0; i1 < 16; i1 += 1) {
for (int i0 = 0; i0 < 16; i0 += 1) {
fieldData_cxiNewLower1[((18*i1)+i0+19)] = - (
4.2426406871193 * fieldData_bLower3[((18*i1)+i0+19)]
* fieldData_cuLower0[((18*i1)+i0+19)]
* fieldData_invDetBLower0[((18*i1)+i0+19)] );
}
}
Enabling optimizations generates better-performing but less readable code, as shown below.
for (int i1 = 0; i1<16; i1 += 1) {
double* const _fieldData_bLower3_p0 =
&(fieldData_bLower3[(18*i1)]);
/* similarly for _fieldData_cxiNewLower1_p0,
_fieldData_cuLower0_p0, _fieldData_invDetBLower0_p0 */
int _start = 0; int _end = 16;
int _intermediate = std::max({_start,(_end-((_end-_start)%4))});
if (_start<_intermediate) {
__m256d _vec00 = _mm256_set1_pd(4.2426406871193);
for (int i0 = _start; i0<_intermediate; i0 += 4) {
__m256d _vec01 =
_mm256_loadu_pd(&(_fieldData_bLower3_p0[(i0+19)]));
__m256d _vec02 =
_mm256_loadu_pd(&(_fieldData_cuLower0_p0[(i0+19)]));
__m256d _vec03 =
_mm256_loadu_pd(&(_fieldData_invDetBLower0_p0[(i0+19)]));
__m256d _vec04;
_vec04 = _mm256_xor_pd(_mm256_mul_pd(
_mm256_mul_pd(_vec02, _vec03),
_mm256_mul_pd(_vec00, _vec01)),
_mm256_set1_pd(-0.0));
_mm256_storeu_pd(&(_fieldData_cxiNewLower1_p0[(i0+19)]),
_vec04);
}
}
for (int i0 = _intermediate; i0<_end; i0 += 1) {
_fieldData_cxiNewLower1_p0[(i0+19)] = - (
4.2426406871193 * _fieldData_bLower3_p0[(i0+19)]
* _fieldData_cuLower0_p0[(i0+19)]
* _fieldData_invDetBLower0_p0[(i0+19)]);
}
}
Block-structured grids¶
To combine the geometric flexibility of unstructured meshes and the performance benefits of unstructured ones, we use block-structured grids (BSGs) for our ocean domains. [ZGAK2019], [ZGAFKK2022] Next to unmasked BSGs, GHODDESS also supports masked ones to accurately mesh fine-scale geometric features and further exploit performance benefits associated with structured grids.
In this example, we use an unmasked BSG with one MPI partition of 58 blocks with 32 elements each. It can be found in GHODDESS-BSG/grids/BSG_bahamas_b1_f58 after checking out the submodule GHODDESS-BSG and is shown below on the left. On the right, we present the bathymetry and the locations of the four recording stations’ exact coordinates in meters
\((38\,667,49\,333)\), \((56\,098,9\,613)\),
\((41\,263,29\,776)\), and \((59\,594,41\,149)\).
Example setup and parameter choice¶
In our example, the bottom friction was given via a standard quadratic friction law \(\tau_{\textrm{bf}} = C_f | \boldsymbol{u} |\) with coefficient \(C_f = 0.009\), and the (constant) Coriolis parameter was set to \(3.19 \times 10^{-5} s^{-1}\).
The tidal forcing was prescribed at the open boundary consisting of five tidal components (O1, K1, N2, M2, S2) given analytically by
with \(t\) (hours) denoting the time elapsed from the beginning of the simulation, and \(\hat{\xi}\) (meters) the specified tidal elevation at the open sea boundary. A so-called cold start initialization is performed: The fluid domain is assumed to be at rest initially (\(\xi_0 =0\), \({\boldsymbol{q}}_0 = 0\)), and the boundary condition is applied gradually over a period of two days using ramping (a coefficient multiplying the amplitude of the tidal forcing set to zero at the beginning of the simulation and linearly increasing to one at the end of day two).
The three-day simulations is run with a constant time step of 15 seconds for the piecewise constant DG discretization (\(p=0\)). Because of the CFL condition, a smaller time step size will be necessary for higher-order approximations or finder grids.
Now we look at the input parameters for this setup.
The parameters are set the file src/config/bahamas.py in the \(\texttt{Config}\) class which interits from \(\texttt{ConfigMaster}\) and \(\texttt{ConfigMasterReal}\) in this case.
Below we describe the most important parameters to be set for running this example:
Time step size and number and approximation order:
config.analytical_test_case = False # to control error output, use physical right hand side and boundary conditions as specified below
config.default_dt = 15 # time step size in seconds
config.use_time_steps = False # True: use a fixed number of time steps (timeSteps), False: use a maximum time (t_max)
config.t_max = 3 * 24 * 60 * 60
config.order = 0 # local approximation order: 0 constant, 1 linear, 2 quadratic, 3 cubic
Grid and bathymetry: Bight of Abaco with Moore’s Island:
config.minLevel = 2
config.maxLevel = config.minLevel # 2*2**maxLevel is the number of triangular elements per dimension
config.grid_folder = "../GHODDESS-BSG/grids/BSG_bahamas_b1_f58"
config.n_blocks = 1 # number of MPI partitions
config.n_frags_per_block = 58 # number of blocks
Initial and boundary condition:
config.exact_eta = lambda pos, t: sp.S(0)
config.exact_u = lambda pos, t: sp.S(0)
config.exact_v = lambda pos, t: sp.S(0)
# use land and open sea boundary condition as specified in the mesh
config.available_boundary_conditions = [BoundaryCondition.Land, BoundaryCondition.OpenSea]
# parameters for open sea boundaries
config.nope = 1 # number of open sea boundaries forcing segments
config.nbfr = 5 # total number of forcing frequencies on open sea boundaries
# forcing frequencies
config.amig = [sp.S(0.000067597751162),
sp.S(0.000072921165921),
sp.S(0.000137879713787),
sp.S(0.000140518917083),
sp.S(0.000145444119418)]
# nodal factor
config.ff = [sp.S(1), sp.S(1), sp.S(1), sp.S(1), sp.S(1)]
# equilibrium argument in degrees for tidal forcing on open sea boundaries
config.face = [deg_to_rad(0), deg_to_rad(0), deg_to_rad(0), deg_to_rad(0), deg_to_rad(0)]
# amplitude of the harmonic forcing function at the open sea boundaries for frequency
config.emo = [sp.S(0.075),
sp.S(0.095),
sp.S(0.1),
sp.S(0.395),
sp.S(0.06)]
# phase (in degrees) of the harmonic forcing function at the open sea boundaries for frequency
config.efa = [deg_to_rad(sp.S(194.806)),
deg_to_rad(sp.S(206.265)),
deg_to_rad(sp.S(340)),
deg_to_rad(sp.S(0)),
deg_to_rad(sp.S(42.9718))]
# center of cpp projection in degrees LONG/LAT
config.slam_0 = deg_to_rad(sp.S(282.5))
config.sfea_0 = deg_to_rad(sp.S(24))
# ramp function option
# 0: no ramp function is used for forcing functions
# 1: a hyperbolic tangent ramp function is specified and applied to the surface elevation specified boundary conditions,
# the tidal potential forcing function as well as the wind and atmospheric pressure forcing functions
config.n_ramp = 1
config.d_ramp = 2 # duration of ramp function (in days)
Right-hand side forcings:
# coriolis option parameter 0: spatially constant, 1: spatially variable
config.ncor = 0
# constant coriolis coefficient
config.cori = sp.S(0.0000319)
# friction option parameter 0: linear, 1: nonlinear
config.noli_bf = 1
# nonlinear bottom friction coefficient, can be dependent on node position
config.cf = sp.S(0.009)
# gravitational acceleration
config.default_g = sp.S(9.81)
Minimum-depth restriction to prevent elements from getting dry:
config.restrict_min_depth = True
config.min_depth_threshold = sp.Rational(1, 4)
config.min_depth_zero_velocities = True # set velocity components to zero in case of min depth restriction
OpenMP parallelization:
config.omp_enabled = True
config.omp_num_threads = 8
vtk and recording station output:
config.print_interval = 3600 # in seconds
config.print_vtk = True
config.stations = [[38666.66, 49333.32],
[56097.79, 9612.94],
[41262.60, 29775.73],
[59594.66, 41149.62]]
config.print_stations_eta = True
config.print_stations_u = True
config.print_stations_v = True
if not config.original_formulation:
config.print_stations_u_tilde = True
config.print_stations_v_tilde = True
Running the code and visualizing results¶
Provided a suitable .platform file in the directory generator, code is run with the command python3 sweMain.py config.bahamas. It automatically generates layer 4 code and subsequently the optimized C++ code with OpenMP as specified. config.bahamas exchanged for other configs as required.
The code is directly executed and outputs the time steps specified in \(\texttt{config.print_interval}\) and the time step size along with a \(L^2\) norm
or an \(L^2\)-error if \(\texttt{config.analytical_test_case = True}\). Furthermore, it outputs a summary of kernel timings at the end.
The generated code can be found in the directory generated/generated. The output files from the simulation are in generated/data.
ParaView can, for example, be used to display the free-surface elevation \(\xi\) or the velocities \(u\) and \(v\) as well as the computational mesh. The video below depicts the approximate solution \(\xi\).
Furthermore, a time series for the respective recording stations can be written and visualized, as shown below for the free-surface elevation \(\xi\) at the first station.
Faghih-Naini, S. Kuckuk, V. Aizinger, D. Zint, R. Grosso and H. Köstler. Quadrature-free discontinuous Galerkin method with code generation features for shallow water equations on automatically generated block-structured meshes. In: Advances in Water Resources 138 (2020) 103552. doi:10.1016/j.advwatres.2020.103552.
Zint, R. Grosso, V. Aizinger and H. Köstler. Generation of block structured grids on complex domains for high performance simulation. In: V.A. Garanzha, L. Kamenski, H. Si. (Eds.), Numerical Geometry, Grid Generation and Scientific Computing. Springer International Publishing, Cham, pp. 87–99. (2020) doi.org/10.1007/978-3-030-23436-2_6.
Zint, R. Grosso, V. Aizinger, S. Faghih-Naini, S. Kuckuk and H. Köstler. Automatic Generation of Load-Balancing-Aware Block-Structured Grids for Complex Ocean Domains. In: T. Robinson, D. Moxey, V. Z. Tomov. (Eds.): Proceedings of the 2022 SIAM International Meshing Roundtable. Virtual Conference, Zenodo, (2022) doi.org/10.5281/zenodo.6562440.