FEniCS Application to a Finite Element Quasi-Geostrophic Model. Linear and non-linear analysis of Munk-like Solutions and Data Assimilation Implementation Through Adjoint Method

TL;DR

This paper develops a finite element model for the quasi-geostrophic equation using FEniCS, compares it with Navier-Stokes, and explores data assimilation via control theory and adjoint methods for large-scale problems.

Contribution

It introduces a finite element approach for quasi-geostrophic equations with FEniCS and demonstrates the application of control theory and adjoint methods for data assimilation.

Findings

Finite element approximation of quasi-geostrophic model successfully implemented.

Comparison shows differences between Navier-Stokes and quasi-geostrophic models.

Data assimilation approach via adjoint method enhances model accuracy.

Abstract

The first aim of this work is to construct a Finite Elements model for quasi-geostrophic equation. This model is analyzed through FEniCS interface. As a second purpose, it shows the potential of the control theory to treat Data Assimilation large scale problems. In section one the finite element approximation of the quasi-geostrophic model is described. In section two some numerical results on different domains are shown. In section three a comparison between Navier Stokes model and quasi-geostrophic one is discussed. Last section is dedicated to a more complex data assimilation/control problem based on quasi-geostrophic equations, to test possible developments and improvements.

Equations

$$\frac{\partial\Delta\psi}{\partial t}+\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}\mathcal{F}(\psi,\Delta\psi)+\frac{\partial\psi}{\partial x}=-\sin(\pi y)+\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\Delta^{2}\psi-\Big{(}\frac{\delta_{S}}{L}\Big{)}\Delta\psi,$$

$$\begin{cases}q=\Delta\psi&\text{ in }\Omega,\\ \displaystyle\frac{\partial q}{\partial t}+\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}\mathcal{F}(\psi,q)+\frac{\partial\psi}{\partial x}-\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\Delta q+\frac{\delta_{S}}{L}q=-\sin(\pi y)&\text{ in }\Omega,\\ q=0&\text{ on }\partial\Omega,\\ \psi=0&\text{ on }\partial\Omega.\end{cases}$$

$$\frac{\partial\Delta\psi}{\partial t}+\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}\mathcal{F}(\psi,\Delta\psi+h(x,y))+\frac{\partial\psi}{\partial x}=-\sin(\pi y)+\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\Delta^{2}\psi-\Big{(}\frac{\delta_{S}}{L}\Big{)}\Delta\psi.$$

$$h(x,y)=e^{-\frac{x}{\frac{{\delta }_I}{L}}}+g(x,y),$$

$$\begin{cases}\displaystyle\int_{\Omega}q\varphi+\int_{\Omega}\nabla\psi\cdot\nabla\varphi=0&\forall\varphi,p\in\mathbb{V},\\ \displaystyle\int_{\Omega}\frac{\partial q}{\partial t}p+\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}\int_{\Omega}(\psi\nabla q\times\mathbf{\hat{k}})\cdot\nabla p+\int_{\Omega}\frac{\partial\psi}{\partial x}p+\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\int_{\Omega}\nabla q\cdot\nabla p+\frac{\delta_{S}}{L}\int_{\Omega}qp=f&\forall\varphi,p\in\mathbb{V},\\ \end{cases}$$

$$(v,w{)}_{\mathbb{V}}={\int }_{\Omega }\nabla v\cdot \nabla w\forall v,w\in \mathbb{V}.$$

$$\begin{cases}\displaystyle\int_{\Omega}q_{\delta}\varphi_{\delta}+\int_{\Omega}\nabla\psi_{\delta}\cdot\nabla\varphi_{\delta}=0&\forall\varphi_{\delta},p_{\delta}\in\mathbb{V_{\delta}},\\ \displaystyle\int_{\Omega}\frac{\partial q_{\delta}}{\partial t}p_{\delta}+\displaystyle\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}\int_{\Omega}(\psi_{\delta}\nabla q_{\delta}\times\mathbf{\hat{k}})\cdot\nabla p_{\delta}\displaystyle+\int_{\Omega}\frac{\partial\psi_{\delta}}{\partial x}p_{\delta}+\\ \qquad\qquad\qquad\qquad\displaystyle+\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\int_{\Omega}\nabla q_{\delta}\cdot\nabla p_{\delta}+\frac{\delta_{S}}{L}\int_{\Omega}q_{\delta}p_{\delta}=f_{\delta}&\forall\varphi_{\delta},p_{\delta}\in\mathbb{V}_{\delta}.\par\end{cases}$$

$$\begin{cases}\displaystyle\frac{\partial\mathbf{u}}{\partial t}-\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\Delta\mathbf{u}+\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}(\mathbf{u}\cdot\nabla)\mathbf{u}+\Big{(}\frac{\delta_{S}}{L}\Big{)}\mathbf{u}+\nabla p=\mathbf{f},\\ \mathrm{div}(\mathbf{u})=0,\end{cases}$$

$$\begin{cases}\displaystyle\int_{\Omega}\frac{\partial\mathbf{u}}{\partial t}\cdot\mathbf{w}+\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}\int_{\Omega}\nabla\mathbf{u}\cdot\nabla\mathbf{w}+\Big{(}\frac{\delta_{I}}{L}\Big{)}^{2}\int_{\Omega}(\mathbf{u}\cdot\nabla)\mathbf{u}\cdot\mathbf{w}+\displaystyle\Big{(}\frac{\delta_{S}}{L}\Big{)}\int_{\Omega}\mathbf{u}\cdot\mathbf{w}-\int_{\Omega}\mathbf{f}\cdot\mathbf{w}-\int_{\Omega}\mathrm{div}(\mathbf{w})p=0\\ \displaystyle\int_{\Omega}\mathrm{div}(\mathbf{u})s=0,\par\end{cases}$$

$$-(1+y)v$$

$$(1+y)u$$

$$-{\int }_{\Omega }(1+y)v{w}_1$$

$${\int }_{\Omega }(1+y)u{w}_2$$

$$\displaystyle\frac{1}{\text{Re}}=\Big{(}\frac{\delta_{I}}{L}\Big{)}^{-2}\Big{(}\frac{\delta_{M}}{L}\Big{)}^{3}=O(10^{-5}).$$

$$J({\psi }_i,f)=\frac{1}{2}{\int }_{\Omega }({\psi }_i-{\bar{\psi }}_i{)}^2+\frac{\alpha }{2}{\int }_{\Omega }{f}^2,$$

$$\underset{({\psi }_i,f)}{\min }J({\psi }_i,f)\text{subjected to }\mathcal{E}({\psi }_i,f)=0,$$

$$\underset{f}{\min }J({\psi }_i(f),f)\text{subjected to }\mathcal{E}({\psi }_i(f),f)=0.$$

Full text

Keywords: Finite Element Method, Quasi-Geostrophic, Navier Stokes, Data Assimilation, Adjoint Method, FEniCS, dolfin-adjoint.

1 Quasi-Geostrophic Model and Numerical

Approximation

In this section we will introduce an unsteady quasi-geostrophic system. It describes a constant density flow moving under Earth’s rotation influence. In particular we will consider the large scale wind-driven ocean circulation. For all physical theory, implications and details we refer to [5, 6, 7].

The equations are numerically treated by a finite element (FE) spatial discretization. Usually, the classical numerical methods used in Geophysical Fluid Dynamics are based on finite differences. This kind of approach is widespread and produces good approximated results, whereas FE methods seem to be discarded, even if they lead to very precise results, most of all where high resolution is needed: flows on coastal boundary layers and on bathymetric constraints. They fit well for post processing, data assimilation and control purposes. We refer to [1, 2] for FE approximation theory.

1.1 Physical Model Strong Formulation

Let us introduce our dynamic system on a bi-dimensional domain $\Omega$. For now $\Omega$ is not specified: in simulations we have used two different geographical domain. They will be introduced in section 2. Let $L=10^{6}$ be the domain length scale. We can express the problem in the following way: to a given choice of a non-linear parameter $\displaystyle\frac{\delta_{I}}{L}$ and dissipative parameters $\displaystyle\frac{\delta_{M}}{L}$ and $\displaystyle\frac{\delta_{S}}{L}$ , find solution $\displaystyle\psi$ that verifies, under no-slip conditions:

[TABLE]

where $\mathcal{F}$ is the Jacobian Determinant operator333Given two differentiable functions $s(x,y)$ and $t(x,y)$, $\displaystyle\mathcal{F}(s,t)=\frac{\partial s}{\partial x}\frac{\partial t}{\partial y}-\frac{\partial s}{\partial y}\frac{\partial s}{\partial x}.$. The equation is non-dimensional, so the parameters $\displaystyle\frac{\delta_{I}}{L},\frac{\delta_{S}}{L},\frac{\delta_{M}}{L}$ have as range $[0,1]$. A classical approach to solve this kind of problem (with a bi-Laplacian operator) is to split the equation, obtaining:

[TABLE]

1.1.1 Bathymetry

To make the model more realistic one can add the influence of bathymetry on the equations. In the quasi-geostrophic case, it is very simple to consider this effect in the equation. Let us suppose that the sea floor can be defined by a smooth function444The value of the function must have the same order of the Rossby number. $h:\Omega\rightarrow\mathbb{R}$. The new state equation to be examined is:

[TABLE]

To describe the Atlantic floor555Atlantic domain experiments are analyzed in subsections 2.2 and 3.2., for example, one could use

[TABLE]

where $g(x,y)$ is a Gaussian function representing the mid-Atlantic Ridge, whereas the exponential function represents steep continental slope.

From the simulations, it turns out that the bathymetry does not influence significantly the results, so in what follows, a flat bottom is assumed.

Now we are ready to treat numerically our problem.

1.2 Physical Model Weak Formulation and Approximation

This subsection describes how to find an approximated weak solution for quasi-geostrophic problem. Variational techniques are not trivial and a deep theoretical mathematical knowledge is needed to fully understand it. We leave details, information and implications to the reader: all the methods that will be mentioned are fully treated in [1, 2].

Suppose $\mathbb{V}=H^{1}_{0}(\Omega)$.666$H^{1}_{0}(\Omega)$ denotes the space of square integrable functions, vanishing on $\partial\Omega$, with square integrable first derivative. After few algebraic treatments, one gets the following formulation: find solutions $\psi$ and $q$ in $\mathbb{V}$ such that:

[TABLE]

where $f=-\sin(\pi y)$.

We endow $\mathbb{V}$ with the inner product:

[TABLE]

Now a discrete approximation $\mathbb{V}_{\delta}\subset\mathbb{V}$ is considered. The discretization $\mathbb{V}_{\delta}$ has been built as a standard FE approximation. The discretized system version is:

[TABLE]

The solutions $q_{\delta}$ and $\psi_{\delta}$ will be good approximations of $q$ and $\psi$ in the inner product norm.

For the computational implementation we used FEniCS interface (https://fenicsproject.org/), exploiting DOLFIN library (https://fenicsproject.org/about/components.html).

2 Numerical Results

This section aims at showing some numerical results in two different domains: the unit square and a North Atlantic Ocean grid777In order to build the North Atlantic mesh, we exploited Freefem++ (http://www.freefem.org/) and gmsh (http://gmsh.info/).. The square is used to demonstrate equation stability. From the results we can also deduce that FE model is comparable to finite differences one. After this analysis, we implemented simulations on a more realistic coastal grid, representing the Northern part of Atlantic Ocean. We underline that time scale of the model is of the order of $O(10^{8})$ seconds.

2.1 Simulation on a Square Domain

As previously introduced, the square $\Omega=[0,1]\times[0,1]$ is considered to be the fluid domain. It represents a wide Ocean surface portion. We ran our model with a time increment $dt=10$ (almost 4 months). Simulation plots at different time steps are shown in the following lines. One can observe that at final time the system reach the steady solution profile. Results have been compared to streamfunction profile presented in [5]. In all the experiments $\displaystyle\frac{\delta_{S}}{L}=0$ is considered. We analysed three cases representing three typical geophysical Ocean behaviours:

• •

linear case, with $\displaystyle\frac{\delta_{I}}{L}=0$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

• •

low non-linear case, with $\displaystyle\frac{\delta_{I}}{L}=\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

• •

highly non-linear case, with $\displaystyle\frac{\delta_{I}}{L}=7\cdot 10^{-2}$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-3}$.

The number of iterations needed to the convergence of the solver is reported. It will be indicated by IN. Time evolution is given by implicit Euler method with null initial condition.

In all the plots cold colors are used to represent the low values of the solution, while hot colors represent high values. Let us focus on the highly non-linear model: final stationary convergence has few variations from [5]. However, our FE solution is physically meaningful and coherent. The results are normalize with respect to maximum values. The simulations can be compared with the Munk-like model presented in [5].

As one can observe, the main difference between linear and non-linear solutions is that the linear system provides a thickening of the current on Western boundary, whereas non-linearity leads the solution to a deviation to North-East.

2.2 Simulation on North Atlantic Ocean

In this subsection simulations are on a more complex and realistic domain. We choose to analyse the streamfunction dynamics on the North Atlantic Ocean. We assume no-slip conditions all around the boundaries and $\displaystyle\frac{\delta_{S}}{L}=0$ is considered. The time step chosen is $dt=10$. We simulated until stationary convergence. Three cases are studied:

• •

linear case, with $\displaystyle\frac{\delta_{I}}{L}=0$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

• •

low non-linear case, with $\displaystyle\frac{\delta_{I}}{L}=\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

• •

highly non-linear case, with $\displaystyle\frac{\delta_{I}}{L}=7\cdot 10^{-2}$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-3}$.

The time evolution and IN index are reported in the following plots. The initial condition for implicit Euler method solver is the null fields. The color map for the streamfunction is as before. As in the previous example, a more pronounced North-East flow from North America to Europe is observed when a higher value of $\displaystyle\frac{\delta_{I}}{L}$ is assumed. The results are normalized with respect to maximum of the solution.

3 Why not Using Navier Stokes Equations?

After this analysis one can move to a more general model for flow dynamics: Navier Stokes equations. Quasi-Geostrophic equations derive directly from them, after some scaling procedures. Quasi-geostrophic system could be seen as an approximation, a simplified model that explain large scale circulation, reliably. By adopting the primitive equations, a more precise description of the process could be, in principle, obtained. But is it convenient? This section describes Navier Stokes general model, its numerical results and how it is instable under high non-linearity hypothesis.

3.1 Strong and Weak Formulation

Let us introduce incompressible Navier Stokes model in a bi-dimensional large scale domain $\Omega$, governed by the following system in strong form:

[TABLE]

where $\mathbf{u}=[u,v]$ and $p$ are flow velocity and pressure, respectively, and $\displaystyle\mathbf{f}=[-\frac{1}{\pi}\cos(\pi y)],0]$888This value is linked to $-\sin(\pi y)$: for details see [5]. is the wind stress. First of all, we are not considering the streamfunction anymore. The current itself is our output. After a little algebra, exploiting no-slip conditions for velocity, the following weak formulation for $\mathbf{u}$ and p is reached: find the two components $[u,v]$ in $\mathbb{V}=H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)$ and $p\in\mathbb{P}=L_{0}^{2}(\Omega)$999$L_{0}^{2}(\Omega)=\left\{p\in L^{2}(\Omega)\;:\;\displaystyle\int_{\Omega}p=0\right\}$, where $L^{2}(\Omega)$ is the space of square integrable functions. such that:

[TABLE]

for all $\mathbf{w}\in\mathbb{V}$ and $s\in\mathbb{P}$. In this new framework, we are not considering Coriolis’ effect. In the strong formulation, it can be recovered, in the $\beta$ plane, adding

[TABLE]

that is, in weak formulation, for all $\mathbf{w}=[w_{1},w_{2}]\in\mathbb{V}$:

[TABLE]

Now, after a finite element approximation, the model is ready to be numerically treated.

3.2 Simulation on a square domain

We considered $\Omega=[0,1]\times[0,1]$, to compare results with quasi-geostrophic model. For the same reason, $\displaystyle\frac{\delta_{S}}{L}=0$ is chosen. A time step $dt=10$ is taken and we ran our simulation exploiting implicit Euler until stationary convergence, from a null initial condition. What we expect is a gyre thickening to the Western boundary for the linear case, whereas a gyre moving to North-East for the non-linear case. We assume no-slip conditions for the velocity field.

• •

linear case, with $\displaystyle\frac{\delta_{I}}{L}=0$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

• •

low non-linear case, with $\displaystyle\frac{\delta_{I}}{L}=\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

The color map is the usual cold-hot low-high. Time steps and IN index are reported. Numerical results are normalized with respect to the maximum of the solution.

Now, let us analyse the highly non-linear configuration with with $\displaystyle\frac{\delta_{I}}{L}=7\cdot 10^{-2}$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-3}$. Navier Stokes equations are very instable. Considering $dt=5$ we have the following result for t = 10.

This plot is not surprising if seen in Reynold number terms. Indeed

[TABLE]

So, for high values of non-linearity, Navier Stokes is an instable model and could not represent large scale wind-ocean circulation in an efficient and reliable way.

3.3 Simulation of the North Atlantic Ocean

One can consider Navier Stokes model to describe circulation in a realistic domain, like the North Atlantic Ocean. The general equations, in principle, are able to better capture flow dynamics, but we have to face a very huge stability issue in the non-linear case, as in the square domain experiment. Using the North Atlantic domain, the algorithm does not converge for $\displaystyle\frac{\delta_{I}}{L}\gg\frac{\delta_{M}}{L}$. This is the reason why highly non-linear case is not reported. We can only show results of:

• •

the linear case, with $\displaystyle\frac{\delta_{I}}{L}=0$ and $\displaystyle\frac{\delta_{M}}{L}=7\cdot 10^{-2}$,

• •

the low non-linear case, with $\displaystyle\frac{\delta_{I}}{L}=\frac{\delta_{M}}{L}=7\cdot 10^{-2}$.

As usual: $\displaystyle\frac{\delta_{S}}{L}=0,dt=10$. Time simulation is based on implicit Euler method, starting from a null fields. No-slip boundary conditions are assumed.

As we expected, in linear case we have a thickening phenomenon on American East coast. In low non-linear case one can see how the current moves towards North-East.

3.4 Some Comments

Let us draw some conclusions on the advantages and disadventages of using the Navier Stokes equations. Indisputably, they are a complete and meaningful tool to describe Ocean dynamics. Although, computationally speaking, they are very difficult to treat. Quasi-geostrophic model can afford better the non-linear issues.

Note that Navier Stokes model and quasi-geostrophic model are comparable in linear and low non-linear cases: the IN index are similar. This underlines how Navier Stokes equations can be an important instrument to study the examples presented. In the case of high non-linearity, stabilized FE methods are needed [3, 12].

4 Data Assimilation

Data Assimilation is a technique that allows to be more reliable in forecasting models. How can we chance model features in order to reach a desired output? What is the forcing term needed in order to have a specific solution? In our case we have built an inverse problem on wind stress101010The idea is based on the scenarios presented in [8].. What is the wind intensity that allows our output to be similar to a predicted scenario? This fits perfectly in Data Assimilation concept: try to incorporate data information in our basic model, changing it, in order to be more reliable.

Data Assimilation is based on optimal control problems, faced through adjoint method. The control theory fundamentals are not taken into account in this work, but they are fully treated in [2, 9, 11]. In [4] adjoint method is specifically applied to ocean circulation models.

4.1 A Specific Example

Let us consider the square domain $\Omega=[0,1]\times[0,1]$ and a quasi-geostrophic description of our dynamics. We are interested to assimilate and control terms on the West boundary layer.

Suppose that we have some dataset $D_{1},\dots,D_{n}$ (representing more powerful simulation data, experimental data, different scenarios…) at some time values $t_{1},\dots,t_{n}$. The interpolation of $D_{1},\dots,D_{n}$ on the considered domain, gives $\bar{\psi}_{1},\dots,\bar{\psi}_{n}$ functions. Let $\psi_{i}$ be the quasi-geostrophic solution at time $t_{i}$ for $i=1,\dots,n$. Our goal is to modify $\psi_{i}$ in order to fit $\bar{\psi}_{i}$. The only way to obtain this result is to change parameters or forcing term. In our case we decided to analyse the wind forcing action needed to reach the desired state solution. For this purpose, we exploited Lagrangian adjoint method. Let us define the following objective functional for all $i$:

[TABLE]

where $\alpha$ is a penalization term. The aim of the problem is to minimize the functional, constrained to the state equations:

[TABLE]

where $\mathcal{E}(\psi_{i},f)=0$ represents our state dynamic.

Supposing that, at least locally, for each $f$ a unique $\psi_{i}$ exists, so we can write our solution as $\psi_{i}(f)$. Under some mathematical hypotheses, the minimization problem is equivalent to find:

[TABLE]

Let $\hat{J}(f):=J(\psi_{i}(f),f)$ be the so called reduced functional. Now we have all the ingredients to solve our issue, thanks to Lagrangian formalism and adjoint method.

4.2 Numerical results

In the following subsection we are going to show some basic results, by applying an adjoint approach to geophysical flow modeling.

First of all, let us clarify the specifics of the experiment. Simulation ran until final time $T=90$. We wanted to enrich our model thanks to given data111111The data have obtained from a finite differences MATLAB simulation. at time $t_{1}=30,\;t_{2}=60,\;t_{3}=90.$ These data were imported in FEniCS interface in order to use dolfin-adjoint (www.dolfin-adjoint.org) library: it is able to derive the discrete adjoint model from a state constrained minimization problem. For the use of this optimization library we refer to [13].

At each assimilation time a $f_{i}^{opt}$ is found and promptly substituted to the old $f$ value. The parametric setting is the following: $\displaystyle\frac{\delta_{S}}{L}=0$ and $\displaystyle\frac{\delta_{I}}{L}=\frac{\delta_{M}}{L}=2\cdot 10^{-2}.$

As usual we normalized the solutions with respect to the maximum value of the solution. As initial wind stress $f=-\sin(\pi y)$ is considered. From the results one can notice that the original solution attempts to reach data information as we expected.

Conclusions and Perspectives

This work underlines that Finite Element approach is comparable to classical finite differences methods. Results and simulations allow to conclude that FE methods could be a valid option to analyse quasi-geostrophic problems. This report has a strong educational purpose and could be an inspiration for other studies.

The first improvement could involve the parametric study of the problem: one can simulate for different values of the parameters and compare the solutions. From the mathematical point of view, parametric problems can be faced trough particular numerical methods, like reduced order methods [10]. These are faster than standard Finite Element approaches. They are very useful when one has to manage a huge amount of parameters. Parametric optimal control problems could be a very powerful instruments to study either geophysical either environmental phenomena and they totally fit in Atmospheric and Oceanographic science.

Acknowledgements

First of all, I want to thank the Instituto Nazionale di Oceanografia e Geofisica Sperimentale (OGS), Trieste, Italy, that allow me to write this report about the results obtained during my period of internship.

In particular, in the development of the ideas presented in this work, I would like to acknowledge the support of Professor Renzo Mosetti, OGS.

For the technical and computational advices, I thank Dr. Francesco Ballarin, SISSA.

I would also show my gratitude to Professor Gianluigi Rozza, SISSA, who introduced me to the control theory applied in marine sciences and incited me in this internship experience.