# Geometry-aligned moving frames by removing spurious divergence in
curvilinear mesh with geometric approximation error

Sehun Chun,
Underwood International College, Yonsei University, South Korea,
<EMAIL_ADDRESS>This research was supported by Brain Research Program
through the National Research Foundation of Korea(NRF) funded by the Ministry
of Science and ICT(NRF-2017M3C7A1049051)

###### Abstract

The vertices of curvilinear elements usually lie on the exact domain. However,
the additional grid points introduced for the higher-order polynomial of
numerical solution often are located slightly out of the domain to generate
geometric approximation error. Considered the main cause of the conservational
error, this error is relatively smaller than discretization error but large
enough to significantly deteriorate a long time integration. This paper
proposes a novel scheme with geometry-aligned moving frames and the removal of
corresponding spurious divergence to reduce the conservational error caused by
geometric approximation error.

##  1 Introduction

A curvilinear mesh has been widely welcomed by various computational
communities because it supports a larger element, particularly in a highly-
curved region. Consequently, it provides a larger time marching step and
better approximation of the domain. Thus, the curvilinear mesh shows better
accuracy and efficiency. However, one critical problem with curvilinear mesh
is the additional grid points beyond the triangular or quadrilateral elements’
vertices. Then are introduced when the order of polynomial for the numerical
solution is higher than one. However, the location of these points does not
always lie exactly on the domain to cause geometric approximation error.
Increasing the polynomial order ($p$) does not reduce this error, and only
decreasing the edge length ($h$) of every element would do.

Several efforts have been made to remove the geometric approximation error in
curvilinear meshes. One of them is a high-order mesh generator NekMesh by
relocating the additional grid points such that they also stay on the curved
domain with the optimal accuracy of integration [8, 9, 13, 14]. The location
of grid points is fixed for each mesh, so is $p$. This is in contrast to the
usual mesh, where no specific value of $p$ is chosen. Without such special
relocation of grid points, it is normally considered that most curvilinear
mesh has a geometric approximation error. The critical effects of the
geometric approximation error for the accuracy of the numerical solution are
considered as follows.

The first effect is the change of metric identity to satisfy the following
equality for a divergence-free vector $\mathbf{F}$ on a curved domain $\Omega$
[6, 12].

$\int_{\Omega}\nabla\cdot\mathbf{F}dx=\sum_{i=1}^{2}(J\boldsymbol{\nu}_{i}F_{i})_{;i}=0,$

where $\boldsymbol{\nu}_{i}$ is the $i$th tangent vector of the curved
coordinate axis $\xi^{i}$, $J$ is the Jacobian, and $;i$ is the covariant
differentiation with respect to axis $\xi^{i}$. In addition to the changes of
$\boldsymbol{\nu}_{i}$ and $\xi^{i}$ which will be described in details as the
next effect, the change of the Jacobian $J$ is caused by geometric
approximation error to deteriorate the metric identity. However, the
mathematical or computational methods to redefine the Jacobian by disregarding
geometric approximation error would be challenging.

The second effect is the change of the tangent vector $\boldsymbol{\nu}_{i}$
and consequently, the surface normal vector $\mathbf{k}$ due to the modified
curved axis $\xi^{i}$. This leads to the immediate projection of the initial
vector, which is defined in the original domain. Thus, the magnitude and
direction of the vector could be altered because the vector lies on a
different tangent plane. More importantly, the covariant derivative of the
vector changes on the mesh with geometric approximation error. For example,
vectors on the protruded grid points from the original domain may generate a
spurious divergence, which should not be generated by vectors in the original
domain. Also, the modified direction of the tangent vector could yield a
different magnitude and direction of differential operators such as divergence
and curl.

This paper proposes a novel scheme mainly to cancel the second effect by
geometric approximation error. The tangent vector of the original domain is
retrieved, and the spurious divergence by the protruded grid points beyond the
original domain is removed. Thus, the proposed scheme achieves better accuracy
and conservational properties, particularly in mass and energy. However, this
scheme is restricted for the case when the first effect of geometric
approximation error remains relatively negligible, i.e., the variation of the
Jacobian $J$ is relatively small such that the metric identity approximately
holds.

The order of this paper is organized as follows. Chapter 2 explains the
changes of divergence by geometric approximation error. Chapter 3 proposes a
mathematical scheme to remove the spurious divergence by geometric
approximation error. In Chapter 4, geometry-aligned moving frames are
explained to construct LOCAL and LOCSPH moving frames, depending on the
alignment of the surface normal vector. Chapter 5 and 6 demonstrate the
numerical tests of static divergence and curl on the sphere along with
corresponding partial differential equations such as conservational laws and
Maxwell’s equations. Chapter 7 displays several numerical tests to show the
accuracy of the proposed scheme, particularly in the conservation of mass and
energy. Discussion is followed in Chapter 8.

Figure 1: Illustration of geometric approximation error and the angle
${\varepsilon}$. Dashed line = surface of the target geometry $\Omega$, solid
line = surface of the domain $\Omega^{\varepsilon}$. The end points
corresponding to the vertices are on the curved domain $\Omega$, but not the
intermediate points.

##  2  Divergence by geometric approximation error

Suppose that $\Omega$ is the curved domain. Suppose that the tessellation of
the domain is obtained as

$\Omega=\cup_{e}\Omega_{e},\hskip
14.22636pt\mbox{where}~{}~{}\Omega_{i}\cap\Omega_{j}=\emptyset$

and all the vertices lie on $\Omega$. If the surface is tessellated with the
Euclidean elements with straight edges, such as regular triangles or
quadrilaterals, the numerical solution with $p=1$ does not require additional
points on the domain. For curvilinear meshes where curved edges are
approximated with a polynomial of $q>2$, the numerical solution of $p>2$
places additional grid points on each element other than the vertices. Some of
the additional grid points of the curvilinear mesh for $p>2$ are not on
$\Omega$ but are slightly misplaced from it. This approximation causes a very
small geometric approximation error which is mostly related to
$p$-approximation. Let us denote the element of the domain with a geometric
approximation error of $\varepsilon$ as $\Omega^{\varepsilon}_{e}$. Also, use
the superscript $\varepsilon$ to represent the corresponding quantity, i.e., a
vector or a differential operator, lies on $\Omega^{\varepsilon}$. Applying
the divergence theorem to an arbitrary vector $\mathbf{v}^{\varepsilon}$ on
$\Omega^{\varepsilon}_{e}$ yields

$\int_{\Omega^{\varepsilon}_{e}}\nabla^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}d\mathbf{x}=\int_{\partial\Omega^{\varepsilon}_{e}}\mathbf{n}^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}ds.$
(2.1)

The flux that enters through the boundaries of the domain is the divergence of
the vector. It is often identified as the source of the flow. One problem with
the above equality is that the geometric error also contributes to the
magnitude of the divergence and consequently to the changes of conservational
properties related to the velocity vector such as mass or energy.

For example, consider the spherical shell of radius $r$ with the metric
$ds^{2}=r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\phi^{2}$. Consider a velocity
vector $\mathbf{v}$ on this sphere such that the velocity is expanded by the
tangent vectors of the two axes of the surface such as
$\mathbf{v}=v_{\theta}\boldsymbol{\theta}+v_{\phi}\boldsymbol{\phi}$. The
exact divergence of this vector on $\Omega$ is given as

$\nabla\cdot\mathbf{v}=\frac{1}{r}\frac{\partial
v_{\theta}}{\partial\theta}+\frac{1}{r\sin\theta}\frac{\partial
v_{\phi}}{\partial\phi}.$

For the domain $\Omega^{\varepsilon}$ where the grid points other than the
vertices do not lie on the exact sphere of radius $r$, the $r$ component of
the vector is not negligible, thus the velocity vector
$\mathbf{v}^{\varepsilon}$ on $\Omega^{\varepsilon}$ is expressed as

$\mathbf{v}^{\varepsilon}=v_{r}^{\varepsilon}\mathbf{r}+({v}_{\theta}-{v}_{\theta}^{\varepsilon})\boldsymbol{\theta}+({v}_{\phi}-{v}_{\phi}^{\varepsilon})\boldsymbol{\phi},$
(2.2)

where $v_{r}^{\varepsilon}$, $v_{r}^{\varepsilon}$, and $v_{r}^{\varepsilon}$
are created by the geometric approximation error. Suppose that the geometric
error is sufficiently small and consequently the vector components created by
the geometric approximation error are correspondingly sufficiently small. The
velocity vector $\mathbf{v}^{\varepsilon}$ is still on the surface, thus if
two independent tangent vectors of $\Omega^{\varepsilon}$ are derived such as
$\mathbf{e}^{\varepsilon}_{1}$ and $\mathbf{e}^{\varepsilon}_{2}$, then
$\mathbf{v}^{\varepsilon}$ is expressed on the surface such as
$\mathbf{v}^{\varepsilon}=v_{1}\mathbf{e}^{\varepsilon}_{1}+v_{2}\mathbf{e}^{\varepsilon}_{2}$.
However, $\mathbf{e}^{\varepsilon}_{1}\cdot\mathbf{r}$ or
$\mathbf{e}^{\varepsilon}_{2}\cdot\mathbf{r}$ is not trivial because of the
geometric approximation error. Then, the divergence of Eq. (2.2) is expressed
as

$\nabla\cdot\mathbf{v}^{\varepsilon}=\frac{1}{r^{2}}\frac{\partial(r^{2}v_{r}^{\varepsilon})}{\partial
r}+\frac{1}{r}\frac{\partial}{\partial\theta}({v}_{\theta}-{v}_{\theta}^{\varepsilon})+\frac{1}{r\sin\theta}\frac{\partial}{\partial\phi}({v}_{\phi}-{v}_{\phi}^{\varepsilon}).$

If the tangent vector of the protruded grid points are aligned parallel to the
curved domain $\Omega$, then it is widely accepted that the value of the grid
points are extended in the direction of the surface normal vector. Thus, no
additional error is generated. If additional grid points causing the geometric
error increases gradually make an angle of $\varepsilon$ with respect to the
exact surface as shown in Fig. 1, then the velocity vector would be projected
on the angled surface to induce the changes of $v_{\theta}^{\varepsilon}$ and
$v_{\phi}^{\varepsilon}$ such as

$\|v_{\theta}^{\varepsilon}\|=\|\mathbf{v}^{\varepsilon}\|\cos\varepsilon\approx\mathcal{O}({\varepsilon}^{2}).$

The same error bound for $v_{\phi}^{\varepsilon}$. However, the change of
$v_{r}^{\varepsilon}$ is of a lower order, equivalent to

$\|v_{r}^{\varepsilon}\|=\|\mathbf{v}^{\varepsilon}\|\sin\varepsilon\approx\mathcal{O}({\varepsilon}).$
(2.3)

For a small $\varepsilon$, the error caused by the geometric approximation
error is thus dominated by $v_{r}^{\varepsilon}$ and its derivative along the
surface normal direction. It means that suppressing the divergence caused by
the geometric approximation error can improve the accuracy significantly.
Therefore, the goal of the scheme is to separate the divergence caused by
geometric approximation error from the exact divergence of the velocity vector
on a curved element with a geometric approximation error
$\Omega^{\varepsilon}_{e}$ as follows.

$\int_{\Omega^{\varepsilon}_{e}}\nabla_{S}\cdot\mathbf{v}^{\varepsilon}d\mathbf{x}+\int_{\Omega^{\varepsilon}_{e}}\nabla_{G}\cdot\mathbf{v}^{\varepsilon}d\mathbf{x}=\int_{\partial\Omega^{\varepsilon}_{e}}\mathbf{n}^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}ds,$
(2.4)

where $\nabla_{G}\cdot$ is the divergence operator caused by the geometric
approximation error and $\nabla_{S}\cdot$ is the difference between the total
divergence and $\nabla_{G}\cdot$. In the previous example, $\nabla_{G}\cdot$
corresponds to the derivatives of $r^{2}v_{r}^{\varepsilon}$ along the surface
normal direction.

##  3 Divergence by Curl Theorem

Consider the following Stokes theorem on the curved element
$\Omega^{\varepsilon}_{e}$ as follows.

$\iint_{\Omega^{\varepsilon}_{e}}(\nabla^{\varepsilon}\times\mathbf{F}^{\varepsilon})\cdot
d\mathbf{S}^{\varepsilon}=\oint_{\partial\Omega^{\varepsilon}_{e}}\mathbf{F}^{\varepsilon}\cdot
d\boldsymbol{\ell}^{\varepsilon}.$ (3.1)

For the curved element $\Omega^{\varepsilon}_{e}$, two orthonormal basis
vectors can be chosen such as $\mathbf{e}^{\varepsilon}_{1}$ and
$\mathbf{e}^{\varepsilon}_{2}$. Then, the surface normal vector can be derived
as
$\mathbf{k}^{\varepsilon}=\mathbf{e}^{\varepsilon}_{1}\times\mathbf{e}^{\varepsilon}_{2}$.
Let us first consider the basis vector $\mathbf{e}^{\varepsilon}_{i}$ lying on
the element $\Omega^{\varepsilon}_{e}$.

By substituting
$\mathbf{F}^{\varepsilon}=\mathbf{k}^{\varepsilon}\times\mathbf{v}^{\varepsilon}$
into Eq. (3.1), we obtain

$\iint_{\Omega^{\varepsilon}_{e}}(\nabla^{\varepsilon}\cdot\mathbf{v}^{\varepsilon})dA-\iint_{\Omega^{\varepsilon}_{e}}\mathcal{G}(\mathbf{k}^{\varepsilon},\mathbf{v}^{\varepsilon})dA=\oint_{\partial\Omega^{\varepsilon}_{e}}{\mathbf{n}}^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}ds,$
(3.2)

where
${\mathbf{n}}^{\varepsilon}=d\boldsymbol{\ell}^{\varepsilon}\times\mathbf{k}^{\varepsilon}$.
The new function $\mathcal{G}(\mathbf{k},\mathbf{v})$ is defined as follows.

$\mathcal{G}(\mathbf{k},\mathbf{v})=\mathbf{k}\cdot(\mathbf{k}\cdot\nabla)\mathbf{v}-\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k},$
(3.3)

where each vector differential operations can be computed as follows.

$\mathbf{C}\cdot(\mathbf{A}\cdot\nabla)\mathbf{B}=\mathbf{C}\cdot[\mathbf{A}\cdot\nabla
B_{x}~{}~{}\mathbf{A}\cdot\nabla B_{y}~{}~{}\mathbf{A}\cdot\nabla B_{z}]^{T}.$

The function $\mathcal{G}(\mathbf{k},\mathbf{v})$ approximately measures the
variation of $\mathbf{v}$ with respect to $\mathbf{k}$ and vice versa. The
first component of $\mathcal{G}(\mathbf{k},\mathbf{v})$ in Eq. (3.3) is always
negligible if $\mathbf{v}$ lies on the tangent plane of $\Omega^{\varepsilon}$
and $\mathbf{k}$ is orthogonal to it. This is because the gradient by the
orthonormal basis $\mathbf{e}_{i}$ lies on the tangent plane and is orthogonal
to the surface normal vector $\mathbf{k}^{\varepsilon}$. For other cases such
that $\mathbf{v}$ and $\mathbf{k}$ are not coincident with $\Omega$, this
component could be not negligible. The second component of
$\mathcal{G}(\mathbf{k},\mathbf{v})$ in Eq. (3.3) represents the variation of
the surface normal vector with respect to the velocity vector
$\mathbf{v}^{\varepsilon}$. This component is not necessarily zero because of
mesh error. This quantity usually converges as the space discretization
increases. Particularly for a high-order method such as spectral h/p method,
$\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k}$ converges exponentially to
zero as will be validated in the following test problems in the next sections.

Eq. (3.2) implies that the flow through the boundaries contributes to the
divergence $\nabla^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}$ and the
geometric source $\mathcal{G}(\mathbf{k},\mathbf{v})$. The divergence operator
$\nabla^{\varepsilon}\cdot$ is also defined on the approximated surface
$\Omega^{\varepsilon}_{e}e$, or more accurately on the tangent plane of
$\Omega^{\varepsilon}$. Thus, the divergence
$\nabla^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}$ in Eq. (3.2) computes all
the divergence by the velocity vector $\mathbf{v}^{\varepsilon}$ on
$\Omega^{\varepsilon}_{e}$. However, it also includes the additional
divergence by the inaccurate geometric approximation. Because
$\mathcal{G}(\mathbf{k},\mathbf{v})$ partly represents the additional
divergence by geometric approximation error, the difference of these two
divergences is closer to the net flux of the flow, i.e.,

$\mbox{Net
flux}\approx\nabla^{\varepsilon}\cdot\mathbf{v}^{\varepsilon}-\mathcal{G}(\mathbf{k},\mathbf{v}).$

Thus, $\mathcal{G}(\mathbf{k},\mathbf{v})$ is called spurious divergence by
geometric approximation error.

In the previous example of Eq. (2.2), the velocity vector
$\mathbf{v}^{\varepsilon}$ on $\Omega^{\varepsilon}_{e}$ has a very small
component $v^{\varepsilon}_{r}$ along the surface normal vector $\mathbf{k}$
of the curved domain $\Omega$. Moreover, the gradient defined by the two
orthonormal bases $\mathbf{e}^{\varepsilon}_{i}$ also has a component along
the direction of $\mathbf{k}$. Thus, the surface normal vector of
$\Omega^{\varepsilon}_{e}$, or $\mathbf{k}^{\varepsilon}$, has a small
component orthogonal to $\mathbf{k}$. Therefore, it is possible that
$(\mathbf{v}^{\varepsilon}\cdot\nabla)\mathbf{k}^{\varepsilon}$ is not
trivial, particularly when the space discretization is relatively coarse.
However, it converges to zero as the space discretization increases. The non-
negligible spurious geometric source $\mathcal{G}(\mathbf{k},\mathbf{v})$
independent of $p$ would emerge, particular when the vector is not on the
tangent plane of $\Omega^{\varepsilon}_{e}$.

##  4 Geometrically-aligned moving frames

On the element with a geometric approximation error
$\Omega^{\varepsilon}_{e}$, vectors and the gradients in the numerical scheme
are defined on the plane expanded by the two basis vectors. This plane is
generally oblique to the tangent planes of $\Omega$. Thus, error occurs in the
representation of the velocity vector and its differentiation. However, the
surface normal vector is always orthogonal to the velocity vector and the
tangent plane of $\Omega^{\varepsilon}_{e}$. Thus,
$\mathbf{k}^{\varepsilon}\cdot(\mathbf{k}^{\varepsilon}\cdot\nabla)\mathbf{v}^{\varepsilon}$
is close to the machine precision, though
$\mathbf{k}^{\varepsilon}\cdot(\mathbf{v}^{\varepsilon}\cdot\nabla)\mathbf{k}^{\varepsilon}$
could be non trivial in a coarse mesh. However,
$\mathbf{k}^{\varepsilon}\cdot(\mathbf{v}^{\varepsilon}\cdot\nabla)\mathbf{k}^{\varepsilon}$
decreases as $p$ increases. Consequently, the magnitude of
$\mathcal{G}(\mathbf{k},\mathbf{v})$ is sufficiently small when refinement
increases, i.e., $h$ decreases or $p$ increases.

Consider the orthonormal basis vectors $\mathbf{e}_{i}$, or moving frames, of
the curved domain $\Omega$. The velocity vector defined on $\Omega$ is fully
expanded by $\mathbf{e}_{i}$ without any projection error. Consider a vector
$\mathbf{F}=\mathbf{k}\times\mathbf{v}$ where $\mathbf{k}$ is the surface
normal vector $\mathbf{k}^{\varepsilon}$ of $\Omega$ and the velocity vector
defined as $\mathbf{v}=v^{1}\mathbf{e}_{1}+v^{1}\mathbf{e}_{2}$. Consequently,
the vector $\mathbf{F}$ also lies on the tangent plane of $\Omega$. However,
the surface normal vector of $\Omega^{\varepsilon}_{e}$ is not orthogonal to
$\mathbf{e}_{i}$ and the velocity vector $\mathbf{v}$. Consequently,
$\mathbf{k}\cdot(\mathbf{k}\cdot\nabla)\mathbf{v}$ is significantly large and
does not converge as $p$ increases. The magnitude of
$\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k}$ remains the same as that
with $\mathbf{e}^{\varepsilon}_{i}$ for $\Omega^{\varepsilon}_{e}$.
Consequently, the magnitude of $\mathcal{G}$ is also significant. The basis
vectors $\mathbf{e}_{i}$ of $\Omega$ for $\Omega^{\varepsilon}_{e}$ are called
geometry-aligned moving frames to stress out that the moving frames
$\mathbf{e}_{i}$ represent the geometry of the original domain $\Omega$.

Geometry-aligned moving frames would be particularly useful when the exact
solution or any parameter of equations for the divergence is given in the
specific coordinate system, such as the spherical coordinate axis. One main
reason that it is not favored is that the numerical scheme with geometry-
aligned moving frames yields more error, particularly in conservational
properties. If moving frames lies on the different tangent plane than one of
the curved element $\Omega^{\varepsilon}_{e}$, the scheme generally produces a
larger conservational error or even more significant discretization error. The
proposed scheme with geometry-aligned moving frames in Eq. (3.2) suggests an
accurate numerical scheme to yield better accuracy and overall accurate
conservational properties, as will be demonstrated in the computational tests
in the next sections.

###  4.1 Two types of LOCAL moving frames

In the following test problems, two types of local orthonormal bases, known as
moving frames, are used. The first is based on the orthonormal bases that lie
on the tangent plane of the surface, which is called as LOCAL moving frames
[4]. In the context of the element-wise mapping from the standard element to
the curved element of the domain, the construction of these bases can be
achieved easily by using the differentiation of the map as follows. Suppose
$s_{1}$ and $s_{2}$ be the two independent axes of the standard element of
unit length $\Omega^{s}$. Suppose a map $\mathbf{c}(s_{1},s_{2})$ from
$\Omega^{s}$ to a curved element $\Omega^{e}$. Note that the map
$\mathbf{c}(s_{1},s_{2})$ represents the coordinate of the curved element
$\Omega^{e}$ and the Jacobian is not constant.

This map corresponds the grid points and the tangent vectors of $\Omega^{s}$
to those of $\Omega^{e}$, respectively. If $\partial/\partial s_{i}$ is
represented as the tangent vector of the domain $\Omega^{s}$, then the map
also sends these tangent vectors onto $\Omega^{e}$ as follows.

$\mathbf{c}\left(\frac{\partial}{\partial
s_{i}}\right)=\frac{\partial\mathbf{c}}{\partial
s_{i}}~{}~{}~{}\in~{}~{}\Omega^{e},$

where the equality holds by the law of vector map. The above equation implies
that the tangent vector of the curved element is obtained by differentiating
the map $\mathbf{c}(s_{1},s_{2})$ with respect to the axis $s_{i},~{}i=1,2$.
The orthonormalization of the derived two tangent vectors, which are not
necessarily orthonormal, creates the local orthonormal basis $\mathbf{e}_{1}$
and $\mathbf{e}_{2}$. The differentiability of these bases is guaranteed by
the sufficient smoothness of the map $\mathbf{c}(s_{1},s_{2})$. The third
basis $\mathbf{e}_{3}$ is obtained by the curl of the two basis vectors. The
obtained LOCAL moving frames are generally discontinuous across the interfaces
of the elements. The additional process could be applied, though it is not
necessary, to construct more Euclidean moving frames by comparing the
divergence of $\mathbf{e}_{1}$ along each edge direction [4].

The LOCAL moving frames coincide with the local shape of the surface and are
thus critically affected by geometric approximation error. This means that
$\mathbf{e}_{3}^{\varepsilon}$ and $\mathbf{k}$ has a very small angle
difference. This angle difference can be used as a quantized value to
represent the geometric approximation error. Geometry-aligned moving frames
could be continuous as the spherical coordinate axis. However, it would be
more advantageous to have discontinuous frames across the interfaces of
elements to avoid the geometric singularities, particularly for a closed
surface such as the sphere.

Then, discontinuous geometry-aligned moving frames can be conveniently
obtained from LOCAL moving frames. The first step is to replace
$\mathbf{e}_{3}$ with $\mathbf{k}/\|\mathbf{k}\|$. The second step is to
project $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$ onto the new tangent space which
has $\mathbf{k}/\|\mathbf{k}\|$ as $\mathbf{e}_{3}$. Another
orthonormalization of the basis vectors creates a set of LOCAL moving frames
without geometric approximation error. When the surface normal vector is the
same as that of the sphere with unit length, we call it as LOCSPH (short for
LOCAL SPHERE) moving frames in this paper.

For the test problems in the next sections, the following spherical mesh is
used. The maximum edge length is $h=0.4$ with 480 curved elements. Each edge
has the maximum order of polynomial of 3. Table 1 shows the failure of the
convergence of mesh error as $p$ increase due to geometric approximation
error.

p | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---
$L_{2}$ | 2.0497e-05 | 3.315e-05 | 3.093e-05 | 3.060e-05 | 3.045e-05
$L_{inf}$ | 7.565e-05 | 1.165e-04 | 1.113e-04 | 1.262e-04 | 1.098e-04
Table 1:  For the exact radius $r_{exact}$ of the domain, mesh error in
$L_{2}$ is defined as
$\sqrt{\sum_{i=1}^{N}{(r_{\textrm{exact}}-{r}_{i})^{2}}/{N}}$ for all the
indices $i$ of the total $N$ grid points. Mesh error in $L_{inf}$ is defined
as $\max_{i}\|r_{\textrm{exact}}-{r}_{i}\|$.

##  5 Divergence Test

Consider two static divergence tests on the sphere of unit length. Let
$x=\sin\theta\cos\phi$, $y=\sin\theta\sin\phi$, and $z=\cos\theta$. Expand the
velocity vector in the spherical coordinate axis such as
$\mathbf{v}=v_{\theta}\boldsymbol{\theta}+v_{\phi}\boldsymbol{\phi}$. The
first test adapts the velocity vector $\mathbf{v}$ defined as follows.

$v_{\theta}=1/\sin\theta,~{}~{}~{}v_{\phi}=1.0.$

The second test uses the well-known Rossby-Haurwitz flow defined such as

$\displaystyle v_{\phi}$
$\displaystyle=\omega\sin\theta+K\sin^{3}\theta(4\cos^{2}\theta-\sin^{2}\theta)\cos
4\phi,$ $\displaystyle v_{\theta}$
$\displaystyle=-4K\sin^{3}\theta\cos\theta\sin 4\phi,$

where $\omega=K=7.848\times 10^{-6}~{}s^{-1}$. The divergence of the both
flows on the sphere is zero, satisfying the metric identities as follows.

$\int_{S^{2}}\nabla\cdot\mathbf{v}^{i}=0,\hskip 14.22636pti=1,~{}2.$

By using Eq. (3.2), a discontinuous Galerkin scheme for the above divergence
with respect to a test function $\varphi$ is expressed as

$\int_{S^{2}}(\nabla\cdot\mathbf{v})\varphi
dA=-\int_{S^{2}}\nabla\varphi\cdot\mathbf{v}dA+\int_{\partial
S^{2}}(\mathbf{n}\cdot\mathbf{v}^{*})\varphi
d\ell-\int_{S^{2}}\mathcal{G}(\mathbf{e^{3}},\mathbf{v})\varphi dA,$ (5.1)

where $*$ represents the corresponding quantity is the numerical flux at the
interfaces. The upwind flux is chosen for its accurate performance.

Fig. 2b displays that the second component
$\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k}$ exponentially converges to
zero as $p$ increases. Therefore, this component in $\mathcal{G}$ is
negligible in a sufficiently fine resolution. However, the magnitude of the
first component of $\mathcal{G}$, or
$\mathbf{k}\cdot(\mathbf{k}\cdot\nabla)\mathbf{v}$, is approximately constant.
This component is significantly large, particularly when moving frames do not
lie on the tangent plane of $\Omega^{\varepsilon}_{e}$ such as the LOCSPH
scheme. The different magnitudes of $\mathcal{G}$ for the LOCAL and LOCSPH
schemes are well presented in Fig. 2a. The positive magnitude implies the
scheme with $\mathcal{G}$ is more accurate.

Table 2 shows the comparison between the $L_{2}$ errors between the scheme of
Eq. (5.1) with the two different local frames, LOCAL and LOCSPH moving frames.
In Table 2, LOCSPH moving frames show better accuracy than LOCAL moving
frames, particularly as $p$ increases. For test 1, the $L_{2}$ error of LOCAL
moving frames is stagnant around $3.94$e$-6$, possibly due to mesh error.
However, the $L_{2}$ error of LOCSPH moving frames converges up to
$2.98$e$-7$. This error is less than 10$\%$ of the error of LOCAL moving
frames. For test 2, the $L_{2}$ error of LOCAL moving frames stagnates around
$1.49$e$-06$, whereas the $L_{2}$ error of LOCSPH moving frames converges up
to $1.47$e$-0.7$.

(a) $\mathbf{k}\cdot(\mathbf{k}\cdot\nabla)\mathbf{v}$

(b) $\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k}$

Figure 2: Convergence of $\mathcal{G}$ for static divergence test. For test 1, $\mathbf{k}\cdot(\mathbf{k}\cdot\nabla)\mathbf{v}$ is approximately constant independent of $p$ as $1.4e$-$16$ and $1.8e$-$3$ for LOCAL and LOCSPH moving frames, respectively. For test 2, it is approximately constant as $6.9e$-$17$ and $7.4e$-$4$ for LOCAL and LOCSPH moving frames, respectively. $\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k}$ exponentially converges to zero as $p$ increases for test 1 and test 2 in LOCAL and LOCSPH moving frames. Test 1 |  | Test 2 |
---|---|---|---
p | LOCAL | LOCSPH |  | LOCAL | LOCSPH
3 | 2.24e-04 | 2.24e-04 |  | 4.33e-05 | 4.32e-05
4 | 2.10e-05 | 2.12e-05 |  | 4.42e-06 | 3.98e-06
5 | 6.98e-06 | 6.26e-06 |  | 1.67e-06 | 4.85e-07
6 | 4.12e-06 | 1.97e-06 |  | 1.48e-06 | 1.62e-07
7 | 3.94e-06 | 3.97e-07 |  | 1.49e-06 | 1.49e-07
8 | 3.94e-06 | 2.98e-07 |  | 1.53e-06 | 1.47e-07
Table 2:  Static divergence test: $L_{2}$ error for test 1 and test 2.

###  5.1 Conservational laws

For a constant velocity vector on the sphere, the following conservational
laws is considered.

$\frac{\partial u}{\partial
t}+\nabla\cdot(u\mathbf{v})=0,~{}~{}~{}\mathbf{v}\in S^{2}.~{}~{}~{}$ (5.2)

Then, the weak formulation of the above conservational law is obtained as

$\int_{S^{2}}\frac{du}{dt}\varphi-\int_{S^{2}}\nabla\varphi\cdot\mathbf{v}udA+\int_{\partial
S^{2}}(\mathbf{n}\cdot\mathbf{v}^{*})u\varphi
d\ell-\int_{S^{2}}\mathcal{G}(\mathbf{e^{3}},u\mathbf{v})\varphi dA=0,$ (5.3)

To test the above conservational laws, a smooth cosine bell is advected around
the sphere. Please refer to the following references for the detailed
description of the test [2, 11, 16]. The original problem is to rotate the
bell just one time. In this test, we rotate it up to ten times to observe the
significant changes of the overall error and conservational error. Time is
marched with the step size of $1.0e$-$04$ by an explicit fourth-order Runge-
Kutta scheme. The cosine bell is rotated at the angle of $\pi/4$ with respect
to the axis of the poles.

Fig. 3 presents error convergences for the cosine bell test problem by Eq.
(5.3). Fig. 3a shows that the overall error of LOCSPH with $\mathcal{G}$
increases as time increases. The differences of the LOCAL and LOCSPH schemes
without $\mathcal{G}$ are presented with the dashed circle ($\Delta$LOCAL) and
dashed triangle($\Delta$LOCSPHnoG), respectively. At $T=10.0$, the difference
is approximately $0.0394\%$ and $0.0198\%$, and at $T=20.0$, the difference is
approximately $0.0568$ and $0.0224$ for the LOCAL and LOCSPH schemes without
$\mathcal{G}$, respectively.

Fig. 3b shows the mass conservational error of LOCSPH with $\mathcal{G}$ and
its differences to other schemes. These Figures validate the better accuracy
of LOCSPHwithG in Eq. (5.3) for conservational laws, in overall error, and
conservational error. At $T$=$10.0$, the difference is approximately $12.0\%$,
and at $T$=$20.0$, the difference is approximately $35.8\%$ for LOCSPH without
$\mathcal{G}$. Conservational error of LOCSPH with $\mathcal{G}$ is less than
a thousand times than that of the LOCAL scheme. Therefore, Fig. 3 demonstrates
the accuracy in the overall error and particularly conservational error.

(a) $L_{2}$ error

(b) Mass conservation error

Figure 3: Conservational laws: error slowly increases as time increases for
$p$=$5$. Every $T$=$2$ completes one turn around the great circle. Solid line
corresponds each error of LOCSPH with $\mathcal{G}$ (LOCSPHwithG). Dashed
lines are error difference to LOCSPH without $\mathcal{G}$, or
$\Delta$LOCSPHnoGA = LOCSPHnoGA - LOCSPHwithG, and difference to that of the
LOCAL scheme, or $\Delta$LOCAL = LOCAL - LOCSPHwithG.

##  6 Curl Test

The curl tests for the curl with respect to the surface normal direction
$(\nabla\times\mathbf{v})\cdot\mathbf{k}$ are similar to the divergence tests.
For the first test, the velocity vector $\mathbf{v}$ is introduced as

$v_{\theta}=1.0,~{}~{}~{}v_{\phi}=1/\sin\theta.$

For the second test, a flipped Rossby-Haurwitz flow between the $\theta$-axis
and $\phi$ axis is used, which is defined such as

$\displaystyle v_{\theta}$
$\displaystyle=\omega\sin\theta+K\sin^{3}\theta(4\cos^{2}\theta-\sin^{2}\theta)\cos
4\phi,$ $\displaystyle v_{\phi}$ $\displaystyle=4K\sin^{3}\theta\cos\theta\sin
4\phi,$

where $\omega$ and $K$ are the same as that of the divergence test. Then, the
both flows are curl-free to have the following equality

$\int_{\Omega}(\nabla\times\mathbf{v})\cdot\mathbf{k}\varphi dA=0.$ (6.1)

Geometric effects on $(\nabla\times\mathbf{v})\cdot\mathbf{k}$ are more
complicated than the divergence of a vector because the magnitude of the curl
only is measured along the surface normal vector. For example, the LOCAL and
LOCSPH have different surface normal direction when their tangent plane could
be different. A discontinuous Galerkin scheme to compute Eq. (6.1) is given as

$\displaystyle\int_{\Omega}(\nabla\times\mathbf{v})\cdot\mathbf{e}_{3}\varphi
dA=-\int_{\Omega}\nabla\varphi\cdot\tilde{\mathbf{v}}^{c}dA+\int_{\partial\Omega}(\mathbf{n}\cdot\tilde{\mathbf{v}}^{c*})\varphi
d\ell$ $\displaystyle+\int
v^{i}\mathbf{e}_{i}\cdot(\nabla\times\mathbf{e}_{3})\varphi
dx+\int_{\Omega}\mathcal{G}(\mathbf{e^{3}},\tilde{\mathbf{v}}^{c})\varphi dA,$
(6.2)

where
$\tilde{\mathbf{v}}^{c}\equiv-(\mathbf{v}\cdot\mathbf{e}_{1})\mathbf{e}_{2}+(\mathbf{v}\cdot\mathbf{e}_{2})\mathbf{e}_{1}$.
Table 3 confirms the convergence similar to the static divergence. The LOCSPH
scheme is implemented with the $\mathcal{G}$ component. The convergence of
$\mathcal{G}$ for the curl tests is similar to that of the static divergence
test, as shown in Fig. 4.

Test 1 |  | Test 2 |
---|---|---|---
p | LOCAL | LOCSPH |  | LOCAL | LOCSPH
3 | 2.24e-04 | 2.24e-04 |  | 3.61e-05 | 3.69e-05
4 | 1.97e-05 | 1.98e-05 |  | 4.12e-06 | 4.24e-06
5 | 7.10e-06 | 7.22e-06 |  | 1.71e-06 | 1.79e-07
6 | 4.202e-06 | 4.36e-06 |  | 1.51e-06 | 1.61e-07
7 | 4.04e-06 | 4.28e-07 |  | 1.51e-06 | 1.62e-07
8 | 4.04e-06 | 4.29e-07 |  | 1.52e-06 | 1.62e-07
Table 3:  Static curl test: $L_{2}$ error for LOCAL and LOCSPH moving frames
versus $p$.  Figure 4: Convergence of $\mathcal{G}$ for the static curl tests.
$\mathcal{G}$ is of a significant magnitude only for LOCSPH moving frames,
independent of $p$. $\mathbf{k}\cdot(\mathbf{k}\cdot\nabla)\mathbf{v}$ is
approximately constant as $1.72e$-$29$ and $5.28e$-$6$ with LOCAL and LOCSPH
moving frames. Similarly, $\mathbf{k}\cdot(\mathbf{v}\cdot\nabla)\mathbf{k}$
exponentially converges as $p$ increases.

###  6.1 Maxwell’s equations

For time-dependent PDEs with the curl operator, Maxwell’s equations are chosen
as follows

$\hat{\varepsilon}\frac{\partial\mathbf{E}}{\partial
t}=\nabla\times\mathbf{H},\hskip
14.22636pt\hat{\mu}\frac{\partial\mathbf{H}}{\partial
t}=-\nabla\times\mathbf{E},$ (6.3)

where $\mathbf{E}$ and $\mathbf{H}$ is the electric field and $h$-field,
respectively. $\hat{\varepsilon}$ and $\hat{\mu}$ is the permittivity and
permeability tensor, respectively. Consider transverse magnetic (TM) mode such
that $\mathbf{H}=H^{1}\mathbf{e}_{1}+H^{2}\mathbf{e}_{2}$ and
$\mathbf{E}=E^{3}\mathbf{e}_{3}$. Transverse electric (TE) mode would yield
the similar results. Then, the weak formulation of Maxwell’s equations in Eq.
(6.3) in the context of discontinuous Galerkin methods is given as follows.
For $i=1,2$, we have

$\displaystyle\int\mu^{i}\frac{\partial H^{i}}{\partial t}\varphi dx$
$\displaystyle=\int\nabla\varphi\cdot\mathbf{E}^{*}_{i}dx-\int\mathbf{E}\cdot(\nabla\times\mathbf{e}_{i})\varphi
dx$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\int_{\partial\Omega}\mathbf{n}\cdot\widetilde{\mathbf{E}}^{*}_{i}\varphi
ds-\int\mathcal{G}(\mathbf{E}^{*}_{i})\varphi dx,$ (6.4)

where $\mathbf{E}=E^{3}\mathbf{e}_{3}$ and
$\mathbf{E}^{*}_{i}={E}^{3}\mathbf{e}_{3i}=E^{3}(\mathbf{e}_{3}\times\mathbf{e}_{i})$.
The tilde sign implies that the corresponding quantity is approximated value
at the boundary, called numerical flux. The permeability tensor is also
decomposed into two frames $\mathbf{e}^{1}$ and $\mathbf{e}^{2}$ such as
$\hat{\mu}=\mu^{1}\mathbf{e}^{1}+\mu^{2}\mathbf{e}^{2}$. The electric field is
updated as follows.

$\displaystyle\int\varepsilon^{3}\frac{\partial E^{3}}{\partial t}\varphi dx$
$\displaystyle=\sum_{i=1}^{2}\left[-\int\nabla\varphi\cdot\mathbf{H}^{*}_{i}dx+\int\mathbf{H}^{i}\cdot(\nabla\times\mathbf{e}_{3})\varphi
dx\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left.\int_{\partial\Omega}\mathbf{n}\cdot\widetilde{\mathbf{H}}^{*}_{i}\varphi
ds+\int\mathcal{G}(\mathbf{H}^{*}_{i})\varphi dx\right],$ (6.5)

where $\mathbf{H}^{i}=H^{i}\mathbf{e}_{i}$ and
$\mathbf{H}^{*}_{i}={H}^{i}\mathbf{e}_{i3}=H^{i}(\mathbf{e}_{i}\times\mathbf{e}_{3})$.
For the test of this numerical scheme, the following test problem is adapted
from ref. [3].

$\frac{\partial\mathbf{H}}{\partial
t}=-\nabla\times\mathbf{E}+\mathbf{F}\sin\omega t,\hskip
14.22636pt\frac{\partial\mathbf{E}}{\partial t}=\nabla\times\mathbf{H},$ (6.6)

where we let $\omega=\sqrt{2}$ and the solution of the above equations is
given as

$\displaystyle\mathbf{H}(\theta,\phi,t)$
$\displaystyle=\mathbf{V}(\theta,\phi)\cos\omega t,\hskip
14.22636pt\mathbf{E}(\theta,\phi,t)=\frac{\zeta(\theta,\phi)}{\omega}\sin\omega
t\mathbf{k},$

where $\mathbf{k}$ is the surface normal vector. The expression of other
variables and the values of parameters are the same as described in ref. [3].
Time step is $dt=0.001$ up to the final time $T=0.5$. A classical Runge Kutta
explicit time step is used for time marching. For the numerical flux, the
efficient upwind flux is used [3].

Table 4 displays the accuracy of the three schemes and confirms again that the
LOCSPH scheme shows the best accuracy. LOCAL moving frames yields the worst
accuracy because the curl operator is computed by basis vectors oblique to the
spherical coordinate axis due to geometric approximation error. Even with a
high $p>3$, the $L_{2}$ error of LOCAL moving frames stagnates around
$1.0e$-$4$. Adapting LOCSPH moving frames improves accuracy, converging up to
$6.86e$-$05$. However, LOCSPH moving frames with $\mathcal{G}$ presents the
best accuracy, which converges up to $4.30e$-$05$.

Fig. 5 presents the initial shock of electromagnetic pulse and energy
conservation error during the extremely low frequency (ELF) propagation on the
sphere. The TM mode is chosen. As shown in Fig. 5a, a shock is initiated at
$(1/\sqrt{2},0,1/\sqrt{2})$ to propagates infinitely between the initial point
and antipode without losing energy. Energy conservation error for various
schemes is shown in Fig. 5b. The LOCAL and LOCSPH schemes provide accurate
conservation of energy. However, LOCSPH with $\mathcal{G}$ demonstrates the
superior performance in the conservation of energy in Maxwell’s equations.
LOCSPH without $\mathcal{G}$ shows an erratic convergence as time marches.
This occurs because moving frames used for the computation of the curl is also
used for the decomposition of vector fields, for example, as shown in Fig.
6.4. Introducing two kinds of moving frames would eliminate this kind of
erratic convergence and will be described in more detail in the next section.

p | LOCAL | LOCSPHnoG | LOCSPHwithG
---|---|---|---
1 | 3.20e-02 | 3.17e-02 | 3.07e-02
2 | 2.61e-03 | 2.62e-03 | 2.61e-03
3 | 5.72e-04 | 4.45e-04 | 4.31e-04
4 | 1.30e-04 | 6.92e-05 | 4.43e-05
5 | 1.53e-04 | 6.86e-05 | 4.30e-05
Table 4: Convergence of error for the problem Eq. (6.6).

(a) Initial shock

(b) Energy error

Figure 5: ELF propagation on the sphere and energy conservation error. TM
mode.

##  7 Application to shallow water equations

An immediate application of the proposed scheme is the shallow water equations
(SWEs) on the sphere given as follows

$\displaystyle\frac{\partial H}{\partial t}+\nabla\cdot(H\mathbf{u})=0,$ (7.1)
$\displaystyle\frac{\partial H\mathbf{u}}{\partial
t}+\nabla\cdot(H\mathbf{u}\mathbf{u})+\frac{g}{2}\nabla
H^{2}=fH(\mathbf{u}\times\mathbf{k})+gH\nabla H_{0},$ (7.2)

where $\eta$ and $\mathbf{u}$ represent the free surface elevation, and the
depth-averaged velocity. $H_{0}$ is the still water depth and $H=H_{0}+\eta$
is the total water depth. $f$ and $g$ means the Coriolis parameter and the
gravitational constant, respectively. Geometry-aligned moving frames for the
SWEs have a different meaning than conservational laws or Maxwell’s equations.
This is because the SWEs are not totally coordinate-independent. A component-
wise expression of Eq. (7.2) for $\mathbf{u}=u\mathbf{e}_{1}+v\mathbf{e}_{2}$
is obtained as follows.

$\displaystyle\frac{\partial H\mathbf{u}}{\partial
t}\cdot\mathbf{e}_{1}+\nabla\cdot(Hu\mathbf{u})+H\mathbf{u}\cdot\nabla\mathbf{e}_{1}\cdot\mathbf{u}+\frac{g}{2}\nabla
H^{2}\cdot\mathbf{e}_{1}=fHu+gH\nabla H_{0}\cdot\mathbf{e}_{1},$ (7.3)
$\displaystyle\frac{\partial H\mathbf{u}}{\partial
t}\cdot\mathbf{e}_{2}+\nabla\cdot(Hv\mathbf{u})+H\mathbf{u}\cdot\nabla\mathbf{e}_{2}\cdot\mathbf{u}+\frac{g}{2}\nabla
H^{2}\cdot\mathbf{e}_{2}=-fHv+gH\nabla H_{0}\cdot\mathbf{e}_{2},$ (7.4)

The above equations assumes the following zero tendency condition.

$\mathbf{u}\cdot\nabla\mathbf{e}_{i}\cdot\mathbf{u}=0,\hskip 28.45274pt1\leq
i\leq 2.$ (7.5)

As moving frames are more Euclidean, the left-hand side of Eq. (7.5) is closer
to zero for $i=1,2$. For the popular Williamson’s tests, the zero tendency of
the spherical coordinate should be satisfied by the spherical coordinate axis.
The tendency is usually not incorporated into the numerical scheme of SWES,
assumed to be zero or approximately negligible compared to discretization
error. Consequently, changing the coordinate axis may change the accuracy,
even on the perfect sphere.

In most numerical schemes [1, 5, 7, 10], the gradient of $H^{2}$, or $\nabla
H^{2}$, is computed in the weak formulation. However, $gH\nabla H_{0}$ is
computed directly because the integration of part for $H\nabla H_{0}$ is
impossible when $H_{0}$ is not a constant. For the given moving frames
$\mathbf{e}_{i}$, the weak formulation of $\nabla H^{2}$ can be computed
similarly as the divergence as follows.

$\int\nabla H^{2}dx=\left(\int(\nabla
H^{2}\cdot\mathbf{e}_{1})\mathbf{e}_{1}\varphi dx\right)+\left(\int(\nabla
H^{2}\cdot\mathbf{e}_{2})\mathbf{e}_{2}\varphi dx\right),$ (7.6)

where the directional derivative can be computed as follows.

$\displaystyle\int\nabla H^{2}\cdot\mathbf{e}_{i}\varphi
dx=-\int\nabla\varphi\cdot\mathbf{e}_{i}H^{2}dx-\int
H^{2}(\nabla\cdot\mathbf{e}_{i})\varphi dx$
$\displaystyle+\int_{\partial\Omega}H^{2}(\mathbf{e}_{i}\cdot\mathbf{n})\varphi
ds+\int\mathcal{G}(H^{2}\mathbf{e}_{i})\varphi dx$ (7.7)

Because of the zero tendency requirement, two sets of moving frames are
introduced for the scheme. The first moving frames, or $\mathbf{e}_{i}$, is
the LOCAL moving frames lying on the tangent plane of the domain. The second
moving frames, or $\mathbf{d}_{i}$, is introduced only to compute the weak
formation of the divergence such as $\nabla\cdot H\mathbf{u}$, $\nabla\cdot
Hu\mathbf{u}$, and $\nabla\cdot Hv\mathbf{u}$. The second moving frames are
aligned along a specific coordinate axis with or without $\mathcal{G}$. Using
Eq (5.1) and (7.7), we obtain the following weak formulation of geometry-
aligned moving frames scheme with $\mathcal{G}$ for the SWE equations in the
context of discontinuous Galerkin method.

$\displaystyle\int\frac{\partial H}{\partial t}\varphi dx-\int
Hu^{d}\mathbf{d}^{1}\cdot\nabla\varphi dx-\int
Hv^{d}\mathbf{d}^{2}\cdot\nabla\varphi dx.$ $\displaystyle\hskip
142.26378pt+\int_{\partial\Omega}\widetilde{H}\tilde{\mathbf{u}}\cdot{\mathbf{n}}\varphi
ds+\int\mathcal{G}(\mathbf{e}_{3},H\mathbf{u})\varphi dx.=0\,,$ (7.8)
$\displaystyle\int\frac{\partial Hu}{\partial t}\varphi dx+\int
H\left(u\frac{\partial\mathbf{e}_{1}}{\partial
t}+v\frac{\partial\mathbf{e}_{2}}{\partial
t}\right)\cdot\mathbf{e}_{1}\varphi\,dx-\int
H{u}(u^{d}\nabla\varphi\cdot\mathbf{d}^{1}+v^{d}\nabla\varphi\cdot\mathbf{d}^{2})dx$
(7.9)
$\displaystyle-\int\frac{gH^{2}}{2}\nabla\varphi\cdot\mathbf{e}_{1}dx+\int_{\partial\Omega}\left[\widetilde{H}\widetilde{u}^{2}+\frac{g}{2}\widetilde{H}^{2}\right](\mathbf{e}_{1}\cdot\mathbf{n})\varphi
dx-\int\frac{gH^{2}}{2}\left(\nabla\cdot\mathbf{e}_{1}\right)\varphi\,dx$
$\displaystyle+\int\mathcal{G}(\mathbf{e}_{3},Hu\mathbf{u}+\frac{g}{2}H^{2}\mathbf{e}_{1})\varphi
dx=\int_{\Omega}f\left(Hv\right)\varphi dx+\int_{\Omega}gH\nabla
H_{0}\cdot\mathbf{e}_{1}\varphi\,dx.$ $\displaystyle\int\frac{\partial
Hv}{\partial t}\varphi dx+\int H\left(u\frac{\partial\mathbf{e}_{1}}{\partial
t}+v\frac{\partial\mathbf{e}_{2}}{\partial
t}\right)\cdot\mathbf{e}_{2}\varphi\,dx-\int
H{v}(u^{d}\nabla\varphi\cdot\mathbf{d}^{1}+v^{d}\nabla\varphi\cdot\mathbf{d}^{2})dx$
(7.10)
$\displaystyle-\int\frac{gH^{2}}{2}\nabla\varphi\cdot\mathbf{e}_{2}dx+\int_{\partial\Omega}\left[\widetilde{H}\widetilde{u}^{2}+\frac{g}{2}\widetilde{H}^{2}\right](\mathbf{e}_{2}\cdot\mathbf{n})\varphi
dx-\int\frac{gH^{2}}{2}\left(\nabla\cdot\mathbf{e}_{2}\right)\varphi\,dx$
$\displaystyle+\int\mathcal{G}(\mathbf{e}_{3},Hv\mathbf{u}+\frac{g}{2}H^{2}\mathbf{e}_{2})\varphi
dx=-\int_{\Omega}f\left(Hu\right)\varphi dx+\int_{\Omega}gH\nabla
H_{0}\cdot\mathbf{e}_{2}\varphi\,dx.$

For the following numerical tests, the Lax-Friedrich flux is used for the
numerical flux, and an explicit fourth-order Runge-Kutta scheme is used for
time marching. The details of the test problems are well described in ref. [3,
16], thus will not be repeated in this paper.

(a) Initial condition

(b) Conservation error

Figure 6: Steady zonal flow: Mass/Energy conservation error. $p=5$. $h=0.4$
with 480 elements. Measured up to $T=5.0$.

###  7.1 Steady Zonal flow

The first test is the steady zonal test with the constant free surface
elevation and velocity vector over time. The angle of the velocity vector with
respect to the polar axis is set as $\pi/4$. The still water depth ($H_{0}$)
is chosen as $1.0$, thus $\nabla H_{0}=0$. We let $\mathbf{d}^{i}$ be LOCSPH
moving frames when $\mathbf{e}_{i}$ is the LOCAL moving frames. Define the
following quantities for convenient comparisons between methods. Moving frames
$\mathbf{e}_{i}$ is chosen LOCAL moving frames. LOCAL$\\_$Mass and
LOCAL$\\_$Energy correspond to the mass and energy conservation error of LOCAL
moving frames for both $\mathbf{e}_{i}$ and $\mathbf{d}^{i}$, respectively.
LOCSPHnoG$\\_$Mass and LOCSPHnoG$\\_$Energy correspond to the mass and energy
conservation error when $\mathbf{d}^{i}$ is LOCSPH without $\mathcal{G}$,
respectively. LOCSPHwithG$\\_$Mass and LOCSPHwithG$\\_$Energy correspond to
the mass and energy conservation error when $\mathbf{d}^{i}$ is LOCSPH with
$\mathcal{G}$, respectively.

$\Delta$LOCAL $\displaystyle=\mbox{LOCAL - LOCSPHwithG},$ $\Delta$LOCSPHnoG
$\displaystyle=\mbox{LOCSPHnoG - LOCSPHwithG}.$

The above quantities are positive if LOCSPHwithG$\\_$Mass and
LOCSPHwithG$\\_$Energy are smaller than the error of LOCAL or LOCSPHwithG.
Note that the plot is in log scale. Thus, the negative quantity does not
appear in the plot.

Figure. 6 displays the increase of the mass and energy conservational error
compared to the initial mass and energy. The LOCSPHwithG scheme confirms the
least conservation error of mass and energy for $h=0.4$ and $p=4$. LOCAL
moving frames show the largest conservational error as its difference is
denoted as the dotted line. The LOCSPHnoG scheme shows approximately 3.5$\%$
and 2.0$\%$ more conservational error than those of the LOCSPHwithG scheme in
mass and energy, respectively. However, this better conservational accuracy of
LOCSPHwithG is not always observed with all $h$ and $p$, especially for a
relatively much finer mesh or relatively much more coarse mesh. It is supposed
that the better accuracy of LOCSPHwithG is observed only when the
conservational error of the divergence component in Eqs. (7.1) and (7.2) is
dominant.

###  7.2 Unsteady Zonal flow

One big difference between the unsteady zonal test and the steady zonal test
is that the still surface elevation($H_{0}$) is not constant as given as
follows.

$H_{0}(\theta)=\frac{133681}{r_{a}g}-\frac{10}{r_{a}g}-\frac{(\Omega_{f}\sin\theta)^{2}}{2g},$

where $r_{a}$ is the radius of the earth, $g$ is the gravitational constant,
and $\Omega_{f}$ is the Earth rotation frequency. Thus, $gH\nabla
H_{0}\cdot\mathbf{e}_{1}$ should be computed directly because it cannot be
computed in a weak form by the integration by part. Consequently, the spurious
diffusion removal scheme cannot be applied to $H\nabla H_{0}$. Consequently,
this component seems to be the main source of conservational error for this
unsteady zonal flow. Fig. 7 confirms the better accuracy of the LOCSPH scheme
with $\mathcal{G}$ compared to LOCAL moving frames. However, it also shows
that the conservational accuracy obtained by $\mathcal{G}$ is quickly
overwhelmed by the conservational error of $gH\nabla
H_{0}\cdot\mathbf{e}_{i}$. The test is run up to $T=0.5$. The LOCSPH scheme
with $\mathcal{G}$ does not show any better accuracy, approximately after
$T=0.35$ in the same mesh as the steady zonal test.

(a) Initial condition

(b) Conservation error

Figure 7: Unsteady zonal flow: Mass/Energy conservation error. $p$=$5$.
$h$=$0.4$ with 480 elements. $L_{2}$ error for LOCAL and LOCSPH with $G$ is
$1.522e$-$05$ and $5.186e$-$06$, respectively.

###  7.3 Rossby-Haurwitz flow

The third test problem is the famous Rossby-Haurwitz flow on the sphere. There
is no exact solution to this problem. However, the periodic distribution of
the free surface elevation ($\eta$) and the velocity vector ($\mathbf{u}$)
along the $\phi$-axis rotate around the axis of the sphere. The conservational
error of mass and energy is considered the most important measure to
demonstrate the accuracy of the implemented numerical schemes. The initial
distribution of $\eta$ and $\mathbf{u}$ is much more complicated than those of
the steady zonal flow, as shown in Fig. 8. However, $H_{0}$ is constant, set
as $1.0$, and $\nabla H_{0}$ is zero. Thus, the LOCSPH with $\mathcal{G}$
scheme can be implemented for better accuracy. Time is run up to 15 days.

Fig. 8 presents the mass and energy conservation error of the LOCSPHwithG
scheme and its differences with the LOCAL and LOCSPHnoG schemes. The general
behavior of the plots is similar to that of the steady zonal flows. The
LOCSPHnoG scheme shows approximately 2.0$\%$ and 1.3$\%$ more conservational
error than those of the LOCSPHwithG in mass and energy, respectively.

(a) Initial condition

(b) Conservation error

Figure 8: Rossy flow: Mass/Energy conservation error. $p=5$. $h=0.4$ with 480
elements.

###  7.4 Perturbed unsteady jet flow

As the last test problem with the constant $H_{0}$, the perturbed unsteady jet
simulation is known to be one of the most challenging tests as the initial
condition is presented in Fig. 9a. The length of the edge should be less than
$h$=$0.1$, and consequently, the number of curved triangular elements are at
least more than $10,000$ for most of the numerical schemes. In such a refined
mesh, the difference between the LOCAL and LOCSPH schemes is negligible. Thus,
the LOCSPHwithG scheme does not always show a better accuracy because the
error of the divergences such as $\nabla\cdot H\mathbf{u}$, $\nabla\cdot
Hu\mathbf{u}$, and $\nabla\cdot Hv\mathbf{u}$ is not dominant in this refined
mesh.

Thus, we test the same problem in a coarser mesh. Note that in this mesh, the
signature pattern of the perturbed unsteady jet flow is not observed due to a
large conservational error. However, it serves the purpose to observe the
accuracy of the LOCSPHwithG scheme. Consider the same sphere as before with
$h=0.4$ with 480 elements. The increasing mass and energy conservation error
of the perturbed jet flow is shown in Fig. 9b. The LOCSPHwithG scheme shows
the best accuracy. The LOCAL moving frames and the LOCSPHnoG scheme shows
approximately a similar conservational error, which is approximately $1\%$
larger than the LOCSPHwithG scheme.

(a) Initial condition

(b) Conservation error

Figure 9: Perturbed unsteady jet: Mass/Energy conservation error. $p=5$.
$h=0.4$ with 480 elements.

###  7.5 On an ellipsoid

Method of moving frames can be applied to any curved surface. Thus, Rossby-
Haurwitz flow was tested again by the SWE equations in Eqs. (7.1) and (7.2) on
an ellipsoid, geometrically similar to the Earth as an ellipsoid defined as
follows [15].

$\left(\frac{x}{a}\right)^{2}+\left(\frac{y}{b}\right)^{2}+\left(\frac{z}{c}\right)^{2}=1,$

where $a/c=b/c=1.003364$ and $c$ is nondimensionalized as the unit length. The
previous tests of SWEs are all derived on the sphere. Thus, we let moving
frames for the divergence LOCSPH, same as the previous tests. Moving frames on
the ellipsoid surface, called LOCAL ellipsoid, are closer to LOCAL moving
frames of this domain. However, LOCAL ellipsoid does not satisfy the time-
invariance conditions of SWEs for steady zonal and Rossby-Haurwitz flow. Thus,
LOCAL ellipsoid produces more conservational error.

On the other hand, this could be considered as a more deformed domain than the
previous cases. In the previous mesh with LOCSPH, the $\mathbf{e}_{3}$ of
LOCAL moving frames at vertices is sufficiently close to $\mathbf{e}_{3}$ of
LOCSPH moving frames. This is because all the vertices lie on the sphere. For
the ellipsoidal domain, the $\mathbf{e}_{3}$ of LOCAL moving frames could be
significantly different from $\mathbf{e}_{3}$ of LOCSPH moving frames. This is
because even at the vertices, the surface normal vector is not aligned along
$\mathbf{e}_{3}$ of LOCSPH moving frames. Thus, this test would show whether
geometry-aligned moving frames can show a better accuracy for the discrete
mesh of a slightly-deformed domain.

Fig. 10 displays the increase of the LOCSPHwithG scheme for the steady zonal
flow and its differences to LOCAL and LOCSPHnoG schemes. For steady-state
flow, as shown in Fig. 10a, the LOCSPHwithG scheme shows the best accuracy
similar to the steady-state flow on the sphere. However, the error accuracy
difference is much larger. The LOCAL scheme has approximately $50\%$ and
$30\%$ more conservational error in mass and energy, respectively. The
LOCSPHnoG scheme has approximately $20\%$ and $10\%$ more conservational error
in mass and energy, respectively.

Fig. 10b displays the increase of mass and energy conservational error for
Rossby-Haurwitz flow. For this complex flow, conservational error from other
components seems to dominate in the beginning. Roughly after $T$=$7$, as more
conservational error accumulates, the LOCSPHwithG scheme seems to successfully
suppress the mass and energy conservational error from the divergence
components to yield better accuracy than LOCAL and LOCSPHnoG schemes.

(a) Steady state

(b) Rossby-Haurwitz flow

Figure 10: On the ellipsoid with the ratio of $a/c=b/c=1.003364$ in
$(x/a)^{2}+(y/b)^{2}+(z/c)^{2}=1$, approximately the same as the Earth, when
$c$ is normalized as the unit length. Mass/Energy conservation error. $p=5$.
$h=0.4$ with 480 elements.

##  8 Discussion

This paper proposes a new scheme to improve the overall accuracy and
particularly conservational error by aligning moving frames along the axis of
the domain and removing spurious divergence. One drawback of the proposed
scheme is that the spurious divergence computation could be expensive,
considering that the obtained accuracy is less than $10\%$ of that of the
original scheme. When the governing equation has multiple components, it is
difficult to know when the divergence error dominates the overall or
conservational error. The proposed scheme also identifies the source of
conservational error as the spurious divergence by geometric approximation
error. This paper specifically implemented the proposed scheme using moving
frames to differentiate and expand vectors on curved surfaces. However, a
similar analysis might apply to other numerical schemes on curvilinear meshes.
The flexibility on the direction of moving frames with spurious divergence
removal, even beyond the tangent plane of the domain, can be adapted in many
future numerical schemes for accuracy, stability, or many other applications.

## References

  * [1] P.-E. Bernard, J.-F. Remacle, R. Comblen, V. Legat, and K. Hillewaert, _High-order discontinuous Galerkin schemes on general 2D manifolds applied to shallow water equations_ , J. Comput. Phys. 228 (2009), no. 17, 6514–6535.
  * [2] Sehun Chun, _Method of moving frames to solve conservation laws on curved surfaces_ , J. Sci. Comput. 53 (2012), no. 2, 268–294.
  * [3] , _Method of moving frames to solve the time-dependent Maxwell’s equations on anisotropic curved surfaces: Applications to invisible cloak and ELF propagation_ , J. Compt. Phys. 340 (2017), 85–104.
  * [4] , _High-order covariant differentiation in application to helmholtz-hodge decomposition on curved surfaces_ , submitted (2020).
  * [5] F. X. Giraldo and T. Warburton, _A nodal triangle-based spectral element method for the shallow water equations on the sphere_ , J. Comput. Phys. 207 (2005), 129–150.
  * [6] David A. Kopriva, _Metric identities and the discontinuous spectral element method on curvilinear meshes_ , J. Sci. Compt. 26 (2006), no. 3, 301–327.
  * [7] M. Läuter, F. X. Giraldo, D. Handorf, and K. Dethloff, _A discontinuous Galerkin method for the shallow water equations in spherical triangular coordinates_ , J. Comput. Phys. 227 (2008), 10226–10242.
  * [8] J. Marcon, M. Turner, D. Moxey, S.J. Sherwin, and J. Peiró, _A variational approach to high-order r-adaptation_ , 26th International Meshing Roundtable, 2017.
  * [9] J. Marcon, M. Turner, J. Peiró, D. Moxey, C. Pollard, H. Bucklow, and M. Gammon, _High-order curvilinear hybrid mesh generation for CFD simulations_ , 2018 AIAA Aerospace Sciences Meeting (Reston, Virginia), no. January, American Institute of Aeronautics and Astronautics, jan 2018.
  * [10] R. D. Nair, S. J. Thomas, and R. D. Loft, _A discontinuous Galerkin global shallow water model_ , Mon. Wea. Rev. 133 (2005), 876–888.
  * [11] Janusz A. Pudykiewicz, _On numerical solution of the shallow water equations with chemical reactions on icosahedral geodesic grid_ , J. Comput. Phys. 230 (2011), 1956–1991.
  * [12] J. F. Thompson, Z. U. Z. Warsi, and C. W. Mastin, _Numerical grid generation_ , Elsevier Science Publishing, 1985.
  * [13] M. Turner, D. Moxey, J. Peiró, M. Gammon, C. Pollard, and H. Bucklow, _A framework for the generation of high-order curvilinear hybrid meshes for CFD simulations_ , Procedia Engineering 203 (2017), 206–218.
  * [14] M. Turner, J. Peiró, and D. Moxey, _Curvilinear mesh generation using a variational framework_ , Computer-Aided Design 103 (2018), 73–91.
  * [15] Wikipedia contributors, _Earth radius — Wikipedia, the free encyclopedia_ , 2020, [Online; accessed 1-October-2020].
  * [16] D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber, _A standard test set for numerical approximations to the shallow water equations in spherical geometry_ , J. Comput. Phys. 102 (1992), no. 1, 211–224.


Second Order parallel tensor on generalized f.pk-space form and hypersurfaces
of generalized f.pk-space form

Punam Gupta and Sanjay Kumar Singh

Abstract. The purpose of the present paper to study a second order symmetric
parallel tensor in generalized f.pk-space form. Second order symmetric
parallel tensor in f.pk-space form is combination of the associated metric
tensor and $1$-forms of structure vector fields. We prove that there does not
exist second order skew-symmetric parallel tensor in f.pk-space form. We also
deduce that there is no parallel hypersurface in a generalized f.pk-space form
but there is semi-parallel hypersurfaces in a generalized f.pk-space form.

MSC2020: 53D10, 53C25, 53C21.

Keywords: f.pk-space form, parallel tensor, parallel hypersurface, semi-
parallel hypersurface.

##  1 Introduction

In 1923, Eisenhart [10] proved that if a positive definite Riemannian manifold
admits a second order parallel symmetric tensor other than a constant multiple
of the metric tensor, then it is reducible. In 1926, Levy [17] proved that a
second order parallel symmetric non-singular (with non-vanishing determinant)
tensor in a real space form is proportional to the metric tensor. As an
improvement of the result of Levy, Sharma [19] proved that any second-order
parallel tensor not necessarily symmetric in a real-space form of dimension
greater than $2$ is proportional to the Riemannian metric. Later on, many
authors [3, 7, 15, 18, 20, 21, 22] studied second-order parallel tensor on
various spaces and obtained the important results. In [13], Gherib and
Belkhelfa studied the second order symmetric parallel tensors on generalized
Sasakian space form and proved that it is proportional to metric tensor.
Recently, Belkhelfa and Mahi [4] proved that second order symmetric parallel
tensors on ${\cal S}$-space form is linear combination of the metric tensor.
and $1$-forms of structure vector fields with constant coefficients. They also
prove the existence of semi-parallel hypersurfaces in ${\cal S}$-space form.

On the other side, Falcitelli and Pastore [11] extend the notion of
generalized Sasakian-space form to the generalized f.pk-space form and
discussed the constancy of the $\varphi$-sectional curvature. They also
interrelate generalized f.pk-space form with generalized Sasakian and
generalized complex space-forms. Also the ${\cal S}$-space form, Sasakian
space form are the particular cases of the generalized f.pk-space form.

The purpose of this paper is to study second order parallel tensors on
generalized f.pk-space form. We prove that second order symmetric parallel
tensors on generalized f.pk-space form is linear combination of the associtaed
positive-definite metric tensor and $1$-forms of structure vector fields. We
prove the non-existence of non-zero second order skew-symmetric parallel
tensor. We also find out that there do not exist a parallel hypersurface of a
generalized f.pk-space forms $M^{2n+s}(F_{1},F_{2},{\cal F})$, tangent to the
structure vector fields. Even there exist semi-parallel hypersurface of a
generalized f.pk-space forms $M^{2n+s}(F_{1},F_{2},{\cal F})$, tangent to the
structure vector fields.

##  2 Preliminaries

An f.pk-manifold [5, 9, 14] is a manifold $M^{2n+s}$ on which is defined an
$f$-structure [23] (also known as $f$-structure with complemented frames or
globally framed $f$-stuctures or $f$-structures with parallelizable kernel
i.e. f.pk-structures), that is a $(1,1)$-tensor field $\varphi$ satisfying

$\varphi^{3}+\varphi=0,$

of rank $2n$ such that the subbundle $\ker\varphi$ is parallelizable. Then
there exists a global frame $\\{\xi_{i}\\},i=1,\ldots,s$ for the subbundle
$\ker\varphi$, with dual $1$-forms $\eta^{i}$, satisfying

$\varphi^{2}=-I+\eta^{i}\otimes\xi_{i},$ (2.1)
$\eta^{i}(\xi_{j})=\delta_{j}^{i},\quad\varphi\,\xi_{i}=0,\quad\eta^{i}\circ\varphi=0.$
(2.2)

An f.pk-structure on a manifold $M^{2n+s}$ is said to be normal if the
Nijenhuis tensor field $N=[\varphi,\varphi]+2d\eta^{i}\otimes\xi_{i}$
vanishes, where $[\varphi,\varphi]$ is the Nijenhuis torsion of $\varphi$.
Consider a Riemannian metric $g$ on $M^{2n+s}$ associated with an f.pk-
structure $(\varphi,\xi_{i},\eta^{i})$ such that

$g\left(\varphi X,\varphi
Y\right)=g\left(X,Y\right)-\sum\limits_{i=1}^{s}\eta^{i}(X)\eta^{i}\left(Y\right),$
(2.3)

for any vector fields $X,Y$ on $\Gamma(TM)$. Then an f.pk-structure is known
as a metric f.pk-structure. A manifold with a metric f.pk-structure is known
as a metric f.pk-manifold.

Let $\Phi$ be the fundamental $2$-form on $M^{2n+s}$ defined by

$\Phi(X,Y)=g(X,\varphi Y)=-g(\varphi X,Y),$

for any vector fields $X,Y$ on $\Gamma(TM)$. A metric f.pk-structure is said
to be a ${\cal K}$-structure [5] if it is normal and the fundamental $2$-form
$\Phi$ is closed, a manifold with a ${\cal K}$-structure is known as ${\cal
K}$-manifold. A ${\cal K}$-structure is said to be a ${\cal S}$-structure if
$d\eta^{i}=\Phi$ for all $i\in\\{1,\ldots,s\\}$, a manifold with an ${\cal
S}$-structure is known as ${\cal S}$-manifold. A ${\cal K}$-structure is said
to be a ${\cal C}$-structure if $d\eta^{i}=0$ for all $i\in\\{1,\ldots,s\\}$,
a manifold with a ${\cal C}$-structure is known as ${\cal C}$-manifold.
manifold. Obviously, if $s=1$, a ${\cal K}$-manifold is a quasi Sasakian
manifold, a ${\cal C}$-manifold is a cosymplectic manifold and a ${\cal
S}$-manifold is a Sasakian manifold.

The Levi-Civita connection $\nabla$ [5, 9] of a metric f.pk-manifold satisfies

$\displaystyle 2g((\nabla_{X}\varphi)Y,Z$ $\displaystyle=$ $\displaystyle
3d\Phi(X,\varphi Y,\varphi Z)-3d\Phi(X,Y,Z)$ $\displaystyle+g(N(Y,Z),\varphi
X)+2d\eta^{j}(\varphi Y,Z)\eta^{j}(X)$ $\displaystyle-2d\eta^{j}(\varphi
Z,Y)\eta^{j}(X)+2d\eta^{j}(\varphi Y,X)\eta^{j}(Z)$
$\displaystyle-2d\eta^{j}(\varphi Z,X)\eta^{j}(Y).$

In particular, for ${\cal S}$-manifolds [5], we have
$\nabla_{X}\xi_{i}=-\varphi X,\quad i=1,\ldots,s$.

A plane section $\Pi$ in $T_{p}M$ spanned by $X$ and $\varphi X$, where $X$ is
a tangent vector orthogonal to structure vector fields, is known as
$\varphi$-section. The sectional curvature of $\varphi$-section is called a
$\varphi$-sectional curvature. A metric f.pk-manifold
$(M^{2n+s},\varphi,\xi_{i},\eta^{i},g)$ has pointwise constant
$\varphi$-sectional curvature if at any point $p\in M^{2n+s},\quad
R_{p}(X,\varphi X,X,\varphi X)$ does not depend on the $\varphi$-section
spanned by $\\{X,\varphi X\\}$.

##  3 Generalized f.pk- space form

###### Definition 3.1

[11] A generalized f.pk- space form $M^{2n+s}(F_{1},F_{2},{\cal F}),$ is a
metric f.pk-manifold $(M^{2n+s},\varphi,\xi_{i},\eta^{i},g)$ which admits
smooth functions $F_{1},F_{2},{\cal F}$ such that its curvature tensor field
satisfies

$\displaystyle R(X,Y)Z$ $\displaystyle=$ $\displaystyle F_{1}\left(g(\varphi
X,\varphi Z)\varphi^{2}Y-g(\varphi Y,\varphi Z)\varphi^{2}X\right)$ (3.1)
$\displaystyle+F_{2}\left(g(Z,\varphi Y)\varphi X-g(Z,\varphi X)\varphi
Y+2g(X,\varphi Y)\varphi Z\right)$
$\displaystyle+\sum\limits_{i,j=1}^{s}F_{ij}\left(\eta^{i}(X)\eta^{j}(Z)\varphi^{2}Y-\eta^{i}(Y)\eta^{j}(Z)\varphi^{2}X\right.$
$\displaystyle\left.+g(\varphi Y,\varphi Z)\eta^{i}(X)\xi_{j}-g(\varphi
X,\varphi Z)\eta^{i}(Y)\xi_{j}\right),$

The $\varphi$-sectional curvature of a generalized f.pk-space form
$M^{2n+s}(F_{1},F_{2},{\cal F})$ is pointwise constant, which is
$c=F_{1}+3F_{2}$.

###### Remark 3.2

For $s\geq 2$, an ${\cal S}$-maniolfd with pointwise constant
$\varphi$-sectional curvature $c$ is known as ${\cal S}$-space form [5, 16].
If the given structure is ${\cal S}$-structure then we obtain a ${\cal
S}$-space form by (3.1) with $F_{1}=\dfrac{c+3s}{4},F_{2}=\dfrac{c-s}{4}$ and
$F_{ij}=1$ for all $i,j\in\\{1,\ldots,s\\}$.

###### Remark 3.3

For $s=1$, we obtain a generalized Sasakian-space form
$M^{2n+1}(f_{1},f_{2},f_{3})$ with $f_{1}=F_{1}$, $f_{2}=F_{2}$ and
$f_{3}=F_{1}-F_{11}$. If the given structure is Sasakian then (3.1) holds with
$F_{11}=1,F_{1}=(c+3)/4,F_{2}=(c-1)/4$ and $f_{3}=F_{1}-F_{11}=(c-1)/4=f_{2}$.
If the given structure is Kenmotsu then (3.1) holds with
$F_{11}=-1,F_{1}=(c-3)/4,F_{2}=(c+1)/4$ and
$f_{3}=F_{1}-F_{11}=(c+1)/4=f_{2}$. If the given structure is cosymplectic
then (3.1) holds with $F_{11}=0,F_{1}=c/4,F_{2}=c/4$ and
$f_{3}=F_{1}-F_{11}=c/4$.

Let $M^{2n+s}$ be a metric f.pk-manifold with its metric tensor $g$ and Levi-
Civita connection $\nabla$. Let $R$ denote the Riemann curvature tensor of
$M^{2n+s}$. If $H$ is a $(0,2)$-tensor which is parallel with respect to
$\nabla$ then we can show easily that

$H(R(X,Y)Z,W)+H(Z,R(X,Y)W)=0.$ (3.2)

###### Theorem 3.4

A second order parallel symmetric tensor in a generalized f.pk- space form
$M^{2n+s}(F_{1},F_{2},{\cal F})$ is a linear combination of the associated
positive definite metric tensor and $1$-forms of structure vector fields if
$s\geq 2$ with assumption that ${\cal F}$ has atleast one non-zero function.

Proof. Using (3.1) in (3.2) and taking $Y=\xi_{k}$ $(k\in\\{1,\ldots s\\})$
and $Z=W$, we get

$H\left(\sum\limits_{j=1}^{s}F_{kj}\left(g(\varphi X,\varphi
Z)\xi_{j}+\eta^{j}(Z)\varphi^{2}X\right),Z\right)=0$ (3.3)

On solving and taking $Z=\xi_{r}(r\in\\{1,\ldots s\\})$, we get

$H(X,\xi_{r})=\sum\limits_{\alpha,\beta=1}^{s}\eta^{\alpha}(X)H(\xi_{\beta},\xi_{r}),\quad
F_{kr}\not=0.$ (3.4)

Using (3.4) in (3.3), we have

$H(X,Z)=g(X,Z)\sum\limits_{\gamma=1}^{s}H(\xi_{j},\xi_{\gamma})-\sum\limits_{\alpha,\beta,\gamma=1}^{s}\eta^{\alpha}(X)\eta^{\beta}(Z)H(\xi_{j},\xi_{\gamma})+\sum\limits_{\alpha,\beta,\gamma=1}^{s}\eta^{\alpha}(X)\eta^{j}(Z)H(\xi_{\beta},\xi_{\gamma}).$
(3.5)

###### Remark 3.5

For $s=1$, the (3.5) reduces to $H(X,Z)=g(X,Z)H(\xi,\xi)$.

###### Remark 3.6

For an ${\cal S}$-maniolfd, $H(\xi_{j},\xi_{\gamma})$ are constant for all
$j,\gamma\in\\{1,\ldots s\\}$. Therefore we can say that a second order
parallel symmetric tensor in a ${\cal S}$-space form is a linear combination
of the associated positive definite metric tensor and $1$-forms of structure
vector fields with constant coefficients if $s\geq 2$ [4].

###### Remark 3.7

For $s=1$ and if the manifold is Sasakian then a second order parallel
symmetric tensor in a Sasakian-space form or generalized Sasakian space form
is proportional to the associated positive definite metric tensor.

As an application of Theorem 3.4, now we consider the Ricci tensor of the
manifold. Since we know that the Ricci tensor $S$ of the manifold is symmetric
$(0,2)$-tensor.

###### Definition 3.8

A non-flat semi-Riemannian manifold $M$ is said to be pseudo-Ricci symmetric
[6] if the Ricci tensor $S$ of type $(0,2)$ of the manifold is non-zero and
satisfies the condition

$(\nabla_{X}S)(Y,Z)=2\alpha(X)S(Y,Z)+\alpha(Y)S(X,Z)+\alpha(Z)S(X,Y),$

for all vector fields $X,Y,Z$, where $\alpha$ is a non-zero $1$-form. If
$\alpha=0$, then the manifold reduces to Ricci symmetric manifold.

###### Definition 3.9

A semi-Riemannian manifold $M$ is said to be Ricci-semisymmetric [8] if its
Ricci tensor $S$ satisfies $R$·$S=0$, that is,

$S(R(X,Y)Z,W)+S(Z,R(X,Y)W)=0$

for all vector fields $X,Y,Z,W$.

###### Definition 3.10

A semi-Riemannian manifold $M$ is said to be Ricci flat if its Ricci tensor
$S$ satisfies $S=0$.

###### Corollary 3.11

A Ricci symmetric generalized f.pk- space form $M^{2n+s}(F_{1},F_{2},{\cal
F})$ is a linear combination of the associated positive definite metric tensor
and $1$-forms of structure vector fields if $s\geq 2$ with assumption that
${\cal F}$ has atleast one non-zero function.

In [1, Theorem 2.2], it was shown that if $M$ is pseudo Ricci symmetric
manifold such that ${\rm div}R=0$, then $M$ is Ricci flat. We know that Ricci
flat condition implies Ricci symmetric but converse need not be. Therefore, we
can state that

###### Corollary 3.12

Let $M^{2n+s}(F_{1},F_{2},{\cal F})$ be a pseudo-Ricci symmetric generalized
f.pk- space form such that ${\cal F}$ has atleast one non-zero function and
$s\geq 2$. Then the linear combination of the associated positive definite
metric tensor and $1$-forms of structure vector fields is zero.

Now, we consider skew-symmetric parallel tensor on generalized f.pk- space
form.

###### Theorem 3.13

There does not exist a non-zero second order parallel skew-symmetric tensor in
a generalized f.pk- space form $M^{2n+s}(F_{1},F_{2},{\cal F})$ with
assumption that ${\cal F}$ has atleast one non-zero function.

Proof. Consider (3.2) with $H$ as a non-zero second order parallel skew-
symmetric tensor, using (3.1) and taking $Y=W=\xi_{k}$ $(k\in\\{1,\ldots
s\\})$ with assumption that ${\cal F}$ has atleast one non-zero function, we
obtain

$H(Z,X)=\eta^{k}(Z)H(X,\xi_{k})+\eta^{k}(X)H(Z,\xi_{k}).$ (3.6)

Let $A$ be the dual $(1,1)$-type tensor which is metrically equivalent to $H$,
that is, $H(X,Y)=g(AX,Y)$. Then (3.6) is equivalent to

$AX=-g(AX,\xi_{k})\xi_{k}+\eta^{k}(X)A\xi_{k}.$ (3.7)

Now, taking the inner product of (3.7) with $\xi_{k}$, we have
$H(X,\xi_{k})=0$. Therefore (3.7) reduces to

$AX=\eta^{k}(X)A\xi_{k}.$ (3.8)

Taking the inner product of (3.8) with $Y$, we have $H(X,Y)=0$.

##  4 Parallel submanifolds of f.pk-space forms

Let $G$ be an immersed hypersurface of $M^{2n+s}(F_{1},F_{2},{\cal F})$. Then
the formulas of Gauss and Weingarten are

$\tilde{\nabla}_{X}Y=\nabla_{X}Y+h(X.Y)N,$ $\tilde{\nabla}_{X}N=-SX,$

where $X$ and $Y$ are tangent vector fields, $N$ a unit normal vector normal
to $G$, $h$ the second fundamental form and $S$ the shape operator of $G$.
Note that $h$ and $S$ are related by $h(X,Y)=g(SX,Y)$. In a hypersurface, the
$(0,4)$ tensor field $\tilde{R}.h$ is defined by

$\tilde{R}.h(X,Y,Z,W)=-h(\tilde{R}(X,Y)Z,W)-h(Z,\tilde{R}(X,Y)W)$ (4.1)

A hypersurface is called parallel [2] if $\tilde{\nabla}h=0$.

###### Theorem 4.1

Let $G$ be an hypersurface of a generalized f.pk-space forms
$M^{2n+s}(F_{1},F_{2},{\cal F})$, tangent to the structure vector fields with
$F_{2}\not=0$, then $G$ is not parallel.

Proof. Let $G$ be a parallel hypersurface of a generalized f.pk-space forms
$M^{2n+s}(F_{1},F_{2},{\cal F})$, tangent to the structure vector fields. Let
$h$ be the second fundamental form of $G$. Consider $N$ be the unit normal of
$G$ in $M$ and $W=-\varphi N$. Since $g(\xi_{i},N)=\eta_{i}(N)=0$ for all $i$,
therefore

$g(W,W)=g(\varphi N,\varphi N)=1,$ $g(N,W)=g(N,\varphi N)=0,$ $\varphi W=N.$

Let $X\in\Gamma(TG)$, then we have

$\varphi X=TX+w(X)N$ (4.2)

where $w$ and $T$ are tensor fields on $G$ of type $(0,1)$ and $(1,1)$,
respectively and $TX$ is also the tangent part of $\varphi X$.

Set $w\not=0$, by (4.2), we have $w(X)=g(X,W)$.

Now, using (3.1), (4.2) and Codazzi equation, we have

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\tilde{\nabla}_{X}h(Y,Z)-\tilde{\nabla}_{Y}h(X,Z)=\left(\tilde{R}(X,Y)Z\right)^{\bot}$
$\displaystyle=$ $\displaystyle F_{2}\left(g(Z,\varphi Y)\varphi X-g(Z,\varphi
X)\varphi Y+2g(X,\varphi Y)\varphi Z\right)^{\bot}$

Since $F_{2}\not=0$, we obtain

$\left(g(Z,\varphi Y)w(X)-g(Z,\varphi X)w(Y)+2g(X,\varphi Y)w(Z)\right)N=0.$

Putting $Z=W$ in above equation, we get

$2g(X,TY)=0,$

which implies that $TY=0$, so $\dim\varphi(T_{p}\left(G\right))=1$ for all
$p\in G$.

Since
$T_{p}\left(M\right)=T_{p}\left(G\right)\oplus\left(T_{p}\left(G\right)\right)^{\bot}$
and $rank\varphi=2n$, therefore we obtain

$2n-1\leq\dim\varphi\left(T_{p}\left(G\right)\right)^{\bot}\leq 2n,$

which is impossible since $n>1$.

###### Remark 4.2

This result is true for $S$-space form [4] and generalized Sasakian-space form
[12].

###### Theorem 4.3

There are semi-parallel hypersurfaces, tangent to the structure vector fields
in a generalized f.pk-space forms $M^{2n+s}(F_{1},F_{2},{\cal F})$ with
assumption that ${\cal F}$ has atleast one non-zero function.

Proof. Let $G$ be not a semi-parallel hypersurfaces, tangent to the structure
vector fields in a generalized f.pk-space forms $M^{2n+s}(F_{1},F_{2},{\cal
F})$. It means that

$\tilde{R}.h(X,Y,Z,W)=-h(\tilde{R}(X,Y)Z,W)-h(Z,\tilde{R}(X,Y)W)\not=0,$

where $h$ is the second fundamental form of $G$, which is symmetric
$(0,2)$-tensor field on $G$. By using the same argument of Theorem 3.4, we get

$H(X,Z)\not=g(X,Z)\sum\limits_{\gamma=1}^{s}H(\xi_{j},\xi_{\gamma})-\sum\limits_{\alpha,\beta,\gamma=1}^{s}\eta^{\alpha}(X)\eta^{\beta}(Z)H(\xi_{j},\xi_{\gamma})+\sum\limits_{\alpha,\beta,\gamma=1}^{s}\eta^{\alpha}(X)\eta^{j}(Z)H(\xi_{\beta},\xi_{\gamma}),$

which is impossible.

###### Remark 4.4

Even for $S$-space form [4], this result is not true.

## References

  * [1] K. Arslan, R. Ezentaş, C. Murathan and C. Özgür, On pseudo Ricci-symmetric manifolds, Balkan J. Geom. Appl. 6 (2001), no. 2, 1-5.
  * [2] C. Asperti, A. Lobos, F. Mercuri, Pseudo-parallel immersions of a space forms, Adv. Geom. 2(2002), no. 1, 57-71.
  * [3] C.L. Bejan and M. Crasmareanu, Parallel second-order tensors on Vaisman manifolds, Int. J. Geom. Methods Mod. Phys. 14 (2017), no. 2, 1750023, 8 pp.
  * [4] M. Belkhelfa and F. Mahi, Second order parallel tensors on ${\cal S}$-manifolds and semi-parallel hypersurfaces of ${\cal S}$-space forms, Ukraïn. Mat. Zh. 71 (2019), no. 10, 1422-1429.
  * [5] D.E. Blair, Geometry of manifolds with structural group $U(n)\times O(s)$, J. Differ. Geom. 4(1970), 155-167.
  * [6] M.C. Chaki, On pseudo-Ricci symmetric manifolds. Bulg. J. Phys. 15 (1988), 525-531.
  * [7] L. Das, Second order parallel tensors on $\alpha$-Sasakian manifold, Acta Math. Acad. Paedagog. Nyházi. (N.S.) 23(2007), no. 1, 65-69.
  * [8] R. Deszcz, On Ricci-pseudosymmetric warped products, Demonstratio Math., 22 (1989), 1053-1065.
  * [9] K. Duggal, S. Ianus and A.M. Pastore, Maps interchanging $f$-structures and their harmonicity, Acta Appl. Math. 67(2001), no. 1, 91-115.
  * [10] L.P. Eisenhart, Symmetric tensors of the second order whose first covariant derivatives are zero, Trans. Amer. Math. Soc. 25 (1923), 297-306.
  * [11] M. Falcitelli and A. M. Pastore, Generalized globally framed $f$-space-forms, Bull. Math. Soc. Sci. Math. Roumanie, 52(2009), no. 3, 291-305.
  * [12] F. Gherib, M. Belkhelfa, Parallel submanifolds of generalized Sasakian space forms, Bull. Transilv. Univ. Braşov Ser. III, 51(2009), no. 2, 185-191.
  * [13] F. Gherib, and M. Belkhelfa, Second order parallel tensors on generalized Sasakian space forms and semi parallel hypersurfaces in Sasakian space forms, Beiträge Algebra Geom. 51 (2010), no. 1, 1-7.
  * [14] S.I. Goldberg and K. Yano, Globally framed $f$-manifolds, Ill. J. Math. 15(1971), 456-474.
  * [15] P. Gupta, Index of pseudo-projectively symmetric semi-Riemannian manifolds, Carpathian Math. Publ. 7(2015), no. 1, 57-65.
  * [16] M. Kobayashi and S. Tsuhiya, Invariant submanifolds of an $f$-manifolds with complemented frames, Kodai Math. Semi. Rep, 24(1972), 430-450.
  * [17] H. Levy, Symmetric tensors of the second order whose covariant derivatives vanish, Annals of Maths. 27(1926), 91-98.
  * [18] A.K. Mondal, U.C. De and C. Özgür, Second order parallel tensors on $\left(\kappa,\mu\right)$-contact metric manifolds, An. Ştiinţ. Univ. ”Ovidius” Constanţa Ser. Mat. 18 (2010), no. 1, 229-238.
  * [19] R. Sharma, Second order parallel tensor in real and complex space forms, Internat. J. Math.Math. Sci. 12(1989), no. 4, 787-790.
  * [20] R. Sharma, Second order parallel tensors on contact manifolds, Algebras Groups Geom. 7(1990), no. 2, 145-152.
  * [21] R. Sharma, Second order parallel tensors on contact manifolds. II, C. R. Math. Rep. Acad. Sci. Canada 13 (1991), no. 6, 259-264.
  * [22] M.M. Tripathi, P. Gupta and J.-S. Kim, Index of quasiconformally symmetric semi-Riemannian manifolds, Int. J. Math. Math. Sci. 2012(2012), 14 pp.
  * [23] K. Yano, On a structure defined by a tensor field $f$ satisfying $f^{3}+f=0$, Tensor N. S. 14(1963), 99-109.

Punam Gupta
Department of Mathematics and Statistics
Dr. Harisingh Gour University
Sagar-470 003, Madhya Pradesh
India
<EMAIL_ADDRESS>

Sanjay Kumar Singh
Department of Mathematics
Indian Institute of Science Education & Research
Bhopal-462 066, Madhya Pradesh
India
<EMAIL_ADDRESS>