\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 24, pp. 1--10.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/24\hfil A stabilized finite element method] {A stabilized finite element method for Stream function vorticity formulation of Navier-Stokes equations} \author[M. Abdelwahed, N. Chorfi, M. Hassine \hfil EJDE-2017/24\hfilneg] {Mohamed Abdelwahed, Nejmeddine Chorfi, Maatoug Hassine} \address{Mohamed Abdelwahed \newline Department of Mathematics, College of Sciences, King Saud University, Riyadh, Saudi Arabia} \email{mabdelwahed@ksu.edu.sa} \address{Nejmeddine Chorfi \newline Department of Mathematics, College of Sciences, King Saud University, Riyadh, Saudi Arabia} \email{nchorfi@ksu.edu.sa} \address{Maatoug Hassine \newline Department of Mathematics, College of Sciences, Monastir University, Tunisia} \email{maatoug\_hassine@yahoo.fr} \dedicatory{Communicated by Vicentiu Radulescu} \thanks{Submitted October 29, 2016. Published January 19, 2017.} \subjclass[2010]{35D30, 76D03, 65N30} \keywords{Stream function; vorticity; stabilized method, \hfill\break\indent high performance computing} \begin{abstract} We the solvability of the two-dimensional stream function-vorticity formulation of the Navier-Stokes equations. We use the time discretization and the method of characteristics order one for solving a quasi-Stokes system that we discretize by a piecewise continuous finite element method. A stabilization technique is used to overcome the loss of optimal error estimate. Finally a parallel numerical algorithm is presented and tested. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}[theorem]{Remark} \newtheorem{corollary}[theorem]{Corollary} \allowdisplaybreaks \section{Introduction} The flow of an incompressible viscous fluid in a two dimensional domain $\Omega$ is characterized by two variables: its velocity $u$ and its pressure $p$. It is described by Navier-Stokes problem \begin{equation}\label{m0} \begin{gathered} \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\nabla \mathbf{u} -\nu \Delta \mathbf{u} +\nabla p = f \quad \text{in } \Omega,\; t>0, \\ \nabla .\mathbf{u} = 0\quad \text{in } \Omega,\; t>0. \end{gathered} \end{equation} This system corresponds to the equation of the conservation of the quantity of movement and the equation of conservation of mass. Here, $\nu$ is the kinematic viscosity of the fluid and $f$ is a given function corresponding to the forces applied to the fluid. In general, we add to this system a boundary condition on the border $\Gamma$ of the domain $\Omega$, like \begin{equation} \label{m7} \mathbf{u}=\mathbf{u}_d \quad \text{on } \Gamma \end{equation} known as the Dirichlet condition. This formulation, called velocity-pressure formulation, can be rewritten by introducing two other scalar functions, called the vorticity (noted by $\omega$) and the stream function (noted by $\psi$) \cite{2,6,113,14,111,CGS,112,222}. The link between these functions is given by the relations: \begin{equation}\label{m1} \omega=\nabla \times \mathbf{u} \quad \text{and}\quad \mathbf{u} =\nabla \times \psi . \end{equation} Then we obtain the following $\psi$-$\omega$ formulation equivalent to \eqref{m0}, \begin{equation}\label{m2} \begin{gathered} \frac{\partial \omega}{\partial t}+ \mathbf{u}\nabla\omega - \nu \Delta \omega = \nabla \times f \quad \text{in } \Omega,\;t>0, \\ \omega +\Delta \psi =0\quad \text{in } \Omega,\; t>0. \end{gathered} \end{equation} The velocity boundary condition \eqref{m7} implies that the function $\psi$ and its normal derivative $\partial_n \psi$ are fixed on the boundary $\Gamma$ and given by $$ \psi= \psi_d \quad\text{and} \quad \frac{\partial \psi}{\partial n} = g \quad \text{on }\Gamma. $$ Such a formulation has two main advantages. The first one is related to the automatic satisfaction of the divergence free condition. The second one concerns the reduction of the number of equations. We propose to solve problem \eqref{m2}. A time discretization of this system using the characteristics method \cite{Russ3,17,333}, leads to study, at each time step $\Delta t$, the following system, called quasi-Stokes problem. \begin{equation} \label{m16} \begin{gathered} \omega + \Delta \psi = 0 \quad \text{in } \Omega \\ -\Delta \omega - \lambda \Delta \psi = F \quad \text{in } \Omega \end{gathered} \end{equation} where $\lambda = \frac{1}{\nu \Delta t}$ and $ F = \lambda \omega ^{P } + \frac{1}{\nu } \nabla \times f $ with $\omega ^{P}(t, x) = \omega (t-\Delta t, x(t-\Delta t))$. The resolution of the system \eqref{m16} by a direct approach leads to the loss of one order in the error estimates. To optimize the behavior of our approach, we adapt the regularization-stabilization technique, introduced in \cite{amara,AD1}, to the Navier-Stokes equations. The main step of the proposed approach is to write the decomposition $\psi-\omega$ in a natural variational framework. Such a decomposition permit us to built a linear piecewise stabilized finite element method having a good behavior and to obtain an optimal error estimate. The numerical resolution of the obtained linear system is performed by the Bi-gradient Conjugate Stabilized method \cite{saa}. After an algorithmic analysis identifying the dependence between the different tasks and data involved, an implementation under MPI (Message Passing Interface) has been done \cite{HPmpi}. We obtain then a parallel stabilized algorithm for the Navier-Stokes problem. We start this paper by presenting respectively the time and spatial discretization of the considered problem. Then we describe the stabilized method and its advantage. The numerical analysis of the obtained system and the validation of the stabilization technique is presented in section 3. In addition we present a parallel algorithm. The performance of the proposed method is illustrated by some numerical results. \section{Discrete problem} We present in this section the time and spatial approximations of problem \eqref{m2}. \subsection{Time discretization} Let $T$ a fixed positive real and $f\in C(0,T,H^{-1}(\Omega))$. We consider the regular partition of $[0,T]$ into $N$ equal subintervals $[t_{i-1},t_i]$, $1\le i \le N$ with $t_0 < t_1 < \ldots < t_N=T$ and $\delta t =t_i-t_{i-1}=\frac{T}{N}$. Let $\psi^{n+1}=\psi(.,t^{n+1})$ and $\omega^{n+1}=\omega(.,t^{n+1})$ being the approximations of $\psi$ and $\omega$ at time $t^{n+1}=(n+1)\delta t$. The time discretization of the problem \eqref{m2} is given by: Find $\psi^{n+1}$ and $\omega^{n+1}$ solutions to \begin{equation}\label{s2n} \begin{gathered} \Big(\frac{\partial \omega}{\partial t}+ \mathbf{u} \nabla \omega\Big)^{n+1} -{\nu} \Delta \omega^{n+1} =\nabla \times f^{n+1} \quad \text{in }\Omega\\ \omega^{n+1} +\Delta \psi^{n+1} =0 \quad \text{in } \Omega\\ \psi^{n+1}= \psi_d \quad \text{on } \Gamma\\ \frac{\partial \psi^{n+1}}{\partial n} = g \quad \text{on } \Gamma \end{gathered} \end{equation} The total derivative of the function omega is approximated using characteristics scheme \cite{Russ3,17} as follows \begin{equation} \big(\frac{\partial }{\partial t}+ \mathbf{u} \nabla \big) \omega (\mathbf x,t^{n+1}) \approx \frac{\omega (\mathbf x,t^{n+1}) -\omega (\mathbf X(\mathbf x,t^{n+1};t^n),t^n)}{\delta t} =\frac{\omega^{n+1}-\omega^n o\mathbf{X}^{n}}{\delta t} \end{equation} where $\mathbf{X}^{n}(\mathbf x)=\mathbf X(\mathbf x,t^{n+1};t^n)$ is the position at time $t^n$ of particle of fluid which is at point $\mathbf x$ at time $t^{n+1}$. System \eqref{s2n} becomes \begin{equation} \label{an9} \begin{gathered} \omega^{n+1}+\Delta \psi^{n+1}=0 \quad \text{in } \Omega \\ \Delta \omega^{n+1} +\frac{1}{\nu \delta t} \Delta \psi^{n+1} = -\frac{1}{\nu}\big(\nabla\times f^{n+1}+\frac{\omega^{n}o\mathbf{X}^{n}}{\delta t} \big) \quad \text{in } \Omega \\ \psi^{n+1}= \psi_d \quad \text{on } \Gamma\\ \frac{\partial \psi^{n+1}}{\partial n} = g\quad \text{on } \Gamma \end{gathered} \end{equation} which is equivalent at each time step to the quasi-Stokes system \begin{equation} \label{an11} \begin{gathered} \omega + \Delta \psi =0 \quad \text{in } \Omega \\ \Delta \omega + \lambda \Delta \psi = F \quad \text{in } \Omega \\ \psi = \psi_d \quad \text{on } \Gamma \\ \frac{\partial \psi}{\partial n} = g \quad \text{on } \Gamma \end{gathered} \end{equation} where $F(\cdot,t^{n+1})=-\frac{1}{\nu} \big(\nabla\times f^{n+1}+\frac{\omega^{n}o\mathbf{X}^{n}}{\delta t}\big)$, and $ \psi_d, g $, $\lambda$, $F$ are given. \subsection{Continuous problem} First, we denote by $ H^{-1}(\Omega)$ the dual space of $$ H_0^1 ( \Omega )= \{\theta \in L^2(\Omega ): \nabla \theta \in L^2 (\Omega ), v_{|_\Gamma} = 0 \} $$ with the associated norm $$ \| v\| _{-1,\Omega} ={\sup_{\varphi \in H_0^1(\Omega)} \frac{\langle v,\varphi \rangle _{-1,1,\Omega}}{ | \varphi | _{1,\Omega}}}, \quad \forall v\in H^{-1}(\Omega) $$ where $\langle\cdot,\cdot\rangle_{-1,1,\Omega}$ is the dual pairing between $ H^{-1}(\Omega)$ and $ H_0^1 ( \Omega )$. We introduce the space $$ H^{-1}(\Delta; \Omega ) = \{ \theta \in L^2(\Omega ):\Delta \theta \in H^{-1}(\Omega )\}, $$ equipped with the norm $$ \| \theta \|_{-1,\Delta ,\Omega } =\big( \| \theta \| _{0,\Omega }^2+\| \Delta \theta \| _{-1,\Omega }^2\big) ^{1/2},\quad \forall \theta \in H^{-1}(\Delta ; \Omega) $$ For simplicity of the analysis, we consider in the following $\psi_d=0$ and $g=0$. A variational formulation of problem \eqref{an11} is given by: Find $(\omega,\psi)\in H^{-1}(\Delta;\Omega)\times H_0^1(\Omega)$ such that \begin{equation} \label{PB1} \begin{gathered} {\int_\Omega\theta \omega\,\mathrm{d}\Omega} + \langle \Delta \theta,\psi\rangle_{-1,1,\Omega} = 0\quad \forall \theta \in H^{-1}(\Delta ; \Omega) \\ \langle\Delta \omega, \eta\rangle _{-1,1,\Omega} - {\lambda \int_\Omega \nabla \eta \nabla \psi \,\mathrm{d}\Omega =\int_\Omega F \eta \,\mathrm{d}\Omega} \quad \forall \eta\in H_0^1(\Omega) \end{gathered} \end{equation} \begin{theorem} \label{thm2.1} For all $ f \in H^{-1}(\Omega)$, problem \eqref{PB1} has a unique solution $(\omega,\psi)\in H^{-1}(\Delta ; \Omega)\times H_0^1(\Omega)$ and there exists a constant $C>0$ such that: \begin{equation} \label{r29} |\psi |_{1,\Omega} + \|\omega \|_{-1,\Delta,\Omega}\leq C \|F \|_{-1,\Omega} \end{equation} Moreover, the solution $(\omega,\psi)$ of \eqref{PB1} is a solution of \eqref{an11}. \end{theorem} For a proof of the above theorem see \cite{3}. \begin{remark}\label{rem1} \rm Using the decomposition \begin{equation}\label{b24} \omega =\omega^0 + \omega^{\ast}, \end{equation} problem \eqref{PB1} can be rewritten as: Find $\omega^0\in H_0^{1}(\Omega)$ such that \begin{equation}\label{b22} \int_{\Omega} {\nabla \omega^0 \nabla \eta }\mathrm{d}\Omega = \langle F,\eta \rangle _{-1,1,\Omega}, \quad \forall \eta \in H_0^{1}(\Omega), \end{equation} and find $\omega^{\ast} \in H^{-1}(\Delta ; \Omega)$, $\psi \in H_0^1(\Omega)$ such that \begin{equation} \label{b23} \begin{gathered} \int_\Omega \theta \omega^{\ast} \,\mathrm{d}\Omega - \int_\Omega \nabla \psi \nabla \theta \,\mathrm{d}\Omega = - \int_\Omega \theta \omega^0 \,\mathrm{d}\Omega \quad \forall \theta \in H^{1}(\Omega), \\ -\langle\Delta \omega^{\ast},\eta \rangle _{-1,1,\Omega} + \lambda \int_\Omega \nabla \psi \nabla \eta \,\mathrm{d}\Omega = 0 \quad \forall \eta \in H_0^{1}(\Omega) \end{gathered} \end{equation} The advantage of this decomposition is to get more regularization on $\omega^{\ast}$, Indeed $-\Delta \omega^{\ast} = \lambda \Delta \psi \in L^2(\Omega)$ while $-\Delta \omega^0 = f \in H^{-1}(\Omega)$ and $-\Delta \omega = f + \lambda \Delta \psi \in H^{-1}(\Omega)$. \end{remark} Next, we shall discretize problem \eqref{b22}-\eqref{b23} using the decomposition defined by \eqref{b24}. This particularity related to $\omega^{\ast}$ will improve the behavior of the discrete method. \subsection{Numerical approximation} Let $(\mathcal{T}_h)_h$ be a regular family of decomposition of $\overline{\Omega}$ in triangles $K$ \cite{13,LZL}. For each triangle $K$, we denote respectively by $h_K$ its diameter and $h={\max_{K\in \mathcal{T}_h}}h_K$. We associate to each decomposition $\mathcal{T}_h$, the set $\mathcal{C}_h$ of the internal edges. Then, for every edge $T$ of $\mathcal{C}_h$, there exists two triangles $K$ and $K'$ of $\mathcal{T}_h$ such that \begin{equation} \label{c3} K\neq K' \quad {and}\quad T=\partial K \cap\partial K'. \end{equation} We define the jump of the normal derivative on each edge $T$ of $\mathcal{C}_h$, by \begin{equation} \label{c4} [ \partial_n v ]_T = \begin{cases} \partial_n^K \,v_{K}+ \partial_n^{K'} \,v_{K'} \text{ p.p. on } T & \text{if } T=\partial K \cap\partial K',\; K, K'\in \mathcal{T}_h \\ 0 \quad \text{on } T & \text{if } T = \partial K \cap \Gamma,\; K \in \mathcal{T}_h , \end{cases} \end{equation} where $v_{K} = v\mid_{K}$ and $ \partial_n^K \,v_{K}$ is the normal derivative of $v_{K}$ on $ \partial K,\, K\in \mathcal{T}_h$, and the discrete scalar product and semi-norm \begin{equation}\label{c9} \begin{gathered} \langle \theta_h \, ,\,\delta_h \rangle_h =\sum_{T \in \mathcal{C}_h}\operatorname{mes}T\int_T [\partial _n\theta_h ] [\partial _n\delta_h ]\,\mathrm{d}T,\quad \forall \theta_h \in X_h, \; \delta_h \in X_h , \\ |\theta_h |_h = \Big(\langle \theta_h ,\,\theta_h \rangle_h \Big)^{1/2} \quad \forall \theta_h \in X_h, \end{gathered} \end{equation} where $\operatorname{mes}T$ denotes the length of the edge of $T$ for $T \in \mathcal{C}_h$. We consider the discrete spaces \begin{gather*} X_h =\{ \theta \in \mathcal{C}^0(\overline{\Omega }): \forall K\in \mathcal{T}_h,\,\theta \mid_K\in \mathbb{P}_1(K)\} \subset H^{-1}(\Delta ; \Omega )\\ X_h^0 = X_h\cap H_0^1(\Omega) \end{gather*} where $\mathbb{P}_1(K)$ is the space of piecewise linear functions defined on $\overline K$. We denote by $\omega_h$, $\psi_h$ and $F_h$ respectively the approximations of $\omega$, $ \psi$ and $F$ in $X_h$. The classical discrete variational formulation of problem \eqref{b22}-\eqref{b23} is written \begin{equation}\label{n3} \begin{gathered} a(\omega_h,\theta_h) + b(\theta_h,\psi_h) = 0 \quad\forall \theta_h \in X_h \\ b(\omega_h,\eta_h) - d(\psi_h,\eta_h) = \int_{\Omega} F\eta_h \,d\Omega \quad\forall \eta_h \in X^o_{h} \\ \psi_h\in X^o_h \quad \omega_h \in X_h, \end{gathered} \end{equation} where $a$, $b$ and $d$ are the following three bilinear forms \begin{gather*} a(\delta_h,\theta_h)=\int_{\Omega}\delta_h \theta_h\,d\Omega ,\quad\forall \delta_h,\, \theta_h \in X_h , \\ b(\theta_h,\varphi_h)=- \int_{\Omega}\nabla \theta_h \nabla \varphi_h \, d\Omega ,\quad\forall \theta_h,\, \varphi_h \in \ X_h ,\\ d(\varphi_h,\eta_h)=\lambda\int_{\Omega}\nabla\varphi_h\nabla\eta_h\,d\Omega, \quad\forall (\varphi_h,\eta_h) \in X_h \times X_h^0. \end{gather*} We remark that the compatibility constant between spaces $X_h$ and $X_h^0$ is independent on $h$ but the associated bilinear form $a(\cdot,\cdot)$ to $X_h$ is elliptic for the $ L^2(\Omega)$ norm. Consequently the coercivity constant depends on $h$. To make up for such inconvenient, we use an approximation method based on a technique of stabilization described in \cite{2,AD1}. It consists in changing this principal form $a(\cdot,\cdot)$ into another form $a_h(\cdot,\cdot)$ by addition a stabilizing term as \begin{equation}\label{19} a_h(\delta_h,\theta_h)=a(\delta_h,\theta_h)+ \beta\,A_h(\delta_h,\theta_h) \quad \forall \delta_h,\;\theta_h \in X_h, \end{equation} with $A_h(\delta_h,\theta_h)=\langle \delta_h,\theta_h\rangle_h$ and $\beta \ge 0$ is a parameter to be chosen. \subsection*{Discrete problem} Using stabilization technique, introduced in \eqref{19}, the discrete variational formulation of \eqref{b22}-\eqref{b23} (which is equivalent to \eqref{PB1}) can be rewritten as \begin{equation}\label{P1h} \begin{gathered} \int_{\Omega} {\nabla \omega^0_h \nabla \eta_h }\mathrm{d}\Omega = \langle f,\eta_h\rangle_{-1,1,\Omega}, \quad \forall \eta_h \in X^{0}_h , \\ \int_{\Omega_h}\omega_h^{\ast}\,\theta_h \,\mathrm{d}{\Omega_h} + \beta\langle \omega_h^{\ast}, \theta_h\rangle_h -\int_{\Omega_h }\nabla \psi_h \nabla \theta_h\,\mathrm{d}\Omega =-\int_{\Omega_h}\omega_h^{0}\,\theta_h \,\mathrm{d}{\Omega_h} \forall \theta_h \in X_h \\ -\int_{\Omega_h} \nabla \omega_h^{\ast}\nabla \eta_h \mathrm{d}{\Omega_h} - \lambda \int_{\Omega_h} \nabla \psi_h\nabla \eta_h\,\mathrm{d}{\Omega_h}=0 \\ \forall \eta_h \in X^{0}_h\; \omega^0_h \in X^0_h,\; \omega^{\ast}_h \in X_h, \psi_h \in X^{0}_h \end{gathered} \end{equation} Setting $ {\omega_h} = \omega^0_h + \omega^{\ast}_h $, the couple $(\psi_h,\omega_h)$ is an approximation of $(\psi,\omega)$. \begin{theorem} \label{thm2.2} The discrete problem \eqref{P1h} has a unique solution. \end{theorem} For a proof of the above theorem see \cite{3}. \subsection*{Error estimates} \begin{theorem}\label{th23} If $\omega^{\ast} \in H^2(\Omega)$ and $\omega^0 \in H^2(\Omega)$, there exist $s \in ]\frac{1}{2},1]$ and $ C >0$, independent of $h$, $\beta$ and $\lambda$, such that \begin{gather} \|\omega^{\ast}-\omega^{\ast}_h \|_{0,\Omega} + \sqrt{\lambda} |\psi - \psi_h |_{1,\Omega} \leq C \big(h^{s+1} | \omega^0 |_{2,\Omega} + e_h \big), \\ \|\omega^0 - \omega^0_h \|_{0,\Omega} \leq , h^{s+1} |\omega^0 |_{2,\Omega}, \\ |\psi - \psi_h |_{1,\Omega} \leq C h |\psi |_{2,\Omega} + C (1 + \sqrt{\beta}) \{ h^{s+1} | \omega^0 |_{2,\Omega} + e_h \} \end{gather} where $$ e_h= \big(\lambda h \sqrt{\beta} +\frac{h}{\sqrt{\beta}} +h \sqrt{\lambda} \big) |\psi |_{2,\Omega} + \big( h^s+\sqrt{\beta} \big) h |\omega^{\ast} |_{2,\Omega} $$ \end{theorem} The proof of the above theorem is an adaptation to the quasi-Stokes problem of the proof presented in \cite{AD1} for the biharmonic problem. We omit it. The parameter $\beta$ is selected in such way that error estimates in theorem \ref{th23} becomes optimal. In the case of $\sqrt{\beta}\,= \frac {1}{\sqrt{\lambda}}$, $$ e_h= h \sqrt{\lambda} |\psi |_{2,\Omega} + \big( h^{s+1} + \frac {h}{\sqrt{\lambda}} \big) |\omega^{\ast} |_{2,\Omega} $$ and we obtain the following error estimates. \begin{corollary} Under hypothesis of theorem \ref{th23}, \begin{gather} \|\omega^{\ast}-\omega^{\ast}_h \|_{0,\Omega} \leq C \{ h \sqrt{\lambda} + h^{s+1} + \frac {h}{\sqrt{\lambda}} \}, \\ |\psi - \psi_h |_{1,\Omega} \leq C \{ h + \frac {h^{s+1}}{\sqrt{\lambda}} + \frac {h}{\lambda} \}, \end{gather} where the constant $C>0$ is independent of $h$ and $\lambda$. \end{corollary} This choice of parameter $\beta$ leads to an optimal error estimates, we get a behavior in $O(h)$. We remark that if $\lambda$ becomes very small, the error estimate on $\psi$ is better than that on $\omega$. \section{Numerical implementation} \subsection{Linear system} Let $n$ the number of mesh nodes, $(\theta_i)_{i=1,n}$ be a base of $ X_h$ and $(\eta_i)_{i=1,n}$ a base of $ X_0^h$. The system \eqref{n3} is written for all $j=1,\dots,n$ \begin{equation} \label{n4} \begin{gathered} \sum_{i=1}^{n}\omega_i\; a(\theta_i,\theta_j) + \beta \sum_{i=1}^{n}\omega_i\; A_h(\theta_i,\theta_j) + \sum_{i=1}^{n}\psi_i \; b(\theta_j,\eta_i)= 0 \\ \sum_{i=1}^{n}\omega_i\; b(\theta_i,\eta_j) - \sum_{i=1}^{n}\psi_i\; d(\eta_i,\eta_j) = \int_\Omega f \eta_j \,d\Omega \end{gathered} \end{equation} which is equivalent to the linear system \begin{equation}\label{n7} M X = N \end{equation} where $$ M=\begin{pmatrix} A + \beta S & \quad B^T\\ B&\quad -D \end{pmatrix}, \quad X=\begin{pmatrix} \omega \\ \psi \end{pmatrix}, \quad N=\begin{pmatrix} 0 \\ F \end{pmatrix} $$ \begin{itemize} \item The matrix $A=(a_{i,j})_{1\le i,j \le n}$ with $a_{ij}=a(\theta_i,\theta_j)$ is computed from the bilinear form $a(\cdot,\cdot)$. \item The matrix $S=(s_{i,j})_{1\le i,j \le n}$ with $s_{ij}=A_h(\theta_i,\theta_j)$ represents the bilinear form $A_h(\cdot,\cdot)$. \item The matrix $B=(b_{i,j})_{1\le i,j \le n}$ with $b_{ij}=b(\eta_i,\eta_j)$. \item $D = -\lambda B$. \item$F$ is associated to the second member with $F_j= \int_\Omega F \eta_j \,d\Gamma $. \end{itemize} For the storage of these matrices, we used the Morse compact method keeping only non zero terms. For the resolution, we implement a parallel Bi-gradient Conjugate Stabilized (BICGSTAB) algorithm \cite{saa}. \subsection{Parallel algorithm} The parallel implementation of the BICGSTAB algorithm was done with fortran and MPI library for communications \cite{HPmpi}. We present in the following some numerical results to show the performance of the parallelized solver. The simulations are performed on a cluster of 8 PC's and for different Meshes (Mesh 1 with $4\times 10^4$ nodes, Mesh 2 with $16\times 10^4$ nodes and Mesh 3 with $64\times 10^4$ nodes). We remark that the number of unknowns is the double of the number of nodes and that from the meshes $1$ to the $2$ and from $2$ to $3$, we resolved a problem with four times more data. \begin{table}[htbp] \caption{Computing time with respect to number of processors} \label{tab1} \begin{center} \begin{tabular}{|c|c|c|c|c|c|} \hline Mesh & 1 proc & 2 proc & 4 proc & 8 proc \\ \hline 1 & 21.05 &11.13 & 6.55 &4.02 \\ \hline 2 & 245.85 &128.04 &70.85 &37.7 \\ \hline 3 &2113.36 & 1061.98 &548.9 &297.65 \\ \hline \end{tabular} \end{center} \end{table} From table \ref{tab1}, we remark that the execution time increases when we raise the number of unknowns and decreases when the number of processors increases. We deduce that more the mesh is finer or more the number of unknown is important, more the gain of time becomes significant. This is proportional to the treated case and is often polluted by the communication time between processors which becomes more significant when the number of processors increases. To evaluate the efficiency of the algorithm, we compute the speed-up which represents the ratio between the execution times for $n$ processors and $1$ processor. The optimal speed-up for $n$ processors is $n$. We present in table \ref{tab2} the obtained speed-up for each mesh. \begin{table}[htbp] \caption{Speed-up}\label{tab2} \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline Mesh & 2 proc & 4 proc & 8 proc \\ \hline %1 & 1.87 &1.01 & 0.63 &0.41 \\ \hline 1 &1.89 & 3.21 &5.25 \\ \hline 2 &1.92 &3.47 &6.52 \\ \hline 3 & 1.99 &3.85 &7.1 \\ \hline \end{tabular} \end{center} \end{table} We remark that the speed-up increases and approaches its optimal when the number of unknowns increases. Indeed, for $2$ processors we obtain a good speed-up for mesh $1$ ($1.89$) which reaches its optimal value ($1.99$) for mesh $3$. On the other hand for $8$ processors we obtain bad speed-up ($5.25$) for mesh $1$ but which becomes better ($7.1$) for mesh $3$. This result is a consequence of the communication time between processors which becomes very important compared to the calculation volume when we increase the number of processors. Therefore, the communication cost degrades the performance of the parallel algorithm. \subsection{Numerical validation} We present in this section a validation of the presented stabilization technique on the quasi-Stokes problem \eqref{an11}. We consider the domain of computation $\Omega=[0,1]\times[0,1]$ and the following analytical stream function solution $$ \psi(x,y) = 3x\sin({\pi x})\cos({\pi y}) \quad 0< x,y < 1. $$ \begin{figure}[htbp] \begin{center} \includegraphics[width=0.4\textwidth, angle=270]{fig1} % err-tot.ps angle=270,width=10cm} \end{center} \caption{Total Error}\label{toterr} \end{figure} \begin{figure}[htbp] \begin{center} \includegraphics[width=0.4\textwidth,angle=270]{fig2} % pente-er.ps ,angle=270 \end{center} \caption{Slope}\label{slope} \end{figure} \begin{figure}[htbp] \begin{center} \includegraphics[width=0.96\textwidth]{fig3a} \\ % diag-om.eps \includegraphics[width=0.32\textwidth]{fig3b} % omega-an.eps \includegraphics[width=0.32\textwidth]{fig3c} % omega-bet0.eps \includegraphics[width=0.32\textwidth]{fig3d} \\ % omega-bet07.eps (a) Exact solution\hfil (b) Without stabilization \hfil (c) With stabilization \end{center} \caption{Vorticity $\omega$.}\label{om-an} \end{figure} \begin{figure}[htbp] \begin{center} \includegraphics[width=0.96\textwidth]{fig4a} \\ % diag-psi.eps \includegraphics[width=0.32\textwidth]{fig4b} % psi-an.eps \includegraphics[width=0.32\textwidth]{fig4c} % psi-bet0.eps \includegraphics[width=0.32\textwidth]{fig4d} \\ % psi-bet07.eps (a) Exact solution \hfil (b) Without stabilization \hfil (c) With stabilization \end{center} \caption{Stream function $\psi$.} \label{ps-an} \end{figure} Figure \ref{toterr} shows the total error $\| \omega _1 -\omega _{1h}\| _{-1,\Delta,\Omega } +| \psi _1 -\psi _{1h}| _{1,\Omega }$ with respect to the stabilization term $\beta$. We remark that the lowest error for this case is obtained for $\beta=0.1$. Figure \ref{slope} shows the slope of the total error according to the parameter $\beta$. We remark that this slope is included in the interval $[1.4,2 ]$ which confirms the obtained error estimate. Figures \ref{om-an} and \ref{ps-an} present a comparison between the exact and computed solutions. Figures \ref{om-an} (a) and \ref{ps-an} (a) show respectively the exact stream function and vorticity isovalues. We give also a comparison between the stabilized method corresponding to $\beta=0.1$ (figures \ref{om-an} (c) and \ref{ps-an} (c)) and the classical finite element method corresponding to $\beta=0$ (figures \ref{om-an} (b) and \ref{ps-an} (b)). One can see the contribution of the stabilization term $\beta$ on the quality of the solution. This contribution is more clear for the vorticity (figure \ref{om-an}(c)). \subsection*{Conclusion} A stabilized finite element method for the Navier-Stokes equations written in stream function-vorticity formulation is presented in this work. In order to optimize the computing time, a parallel implementation for the obtained solver is presented. The proposed approach will be extended in a forthcoming work to a three dimensional Navier-Stokes equations. we denote that the proposed algorithm is general and can be used for many engineering problems related to the Navier-Stokes equations. \subsection*{Acknowledgements} The authors would like to extend their sincere appreciation to the Deanship of Scientific Research at King Saud University for funding this Research group No (RG-1435-026). \begin{thebibliography}{00} \bibitem{amara} M. Amara; \emph{Une m\'ethode optimale de classe C$^{\circ}$ d'approximation du bilaplacien}, Comptes Rendus de l'Academie des Sciences, 319 (1994), 1327-1330. \bibitem{2} M. Amara, C. Bernardi; \emph{Convergence of a finite element discretization of the Navier-Stokes equations in vorticity and stream function formulation}, Mathematical Modelling and Numerical analysis, 33(5) (1999) ,1033-1056. \bibitem{AD1} M. Amara, F. Dabaghi; \emph{An optimal C$^o$ finite element method for the 2D biharmonic problem}, {Numerisch Mathematik}, 90(1) (2001), 19-46. \bibitem{3} F. Brezzi, M. Fortin; \emph{Mixed and hybrid finite element method}, Springer Series in Computational Mathematics 15, Springer Verlag, New York, 1991. \bibitem{6} C. Bernardi, V. Girault, Y. Maday; \emph{Mixed spectral element approximation of the Navier Stokes equations in the stream function and vorticity formulation}, IMA Journal of Numerical Analysis, 12 (1992), 565-608. \bibitem{CGS} P. Colli, G. Gilardi, J. Sprekels; \emph{A boundary control problem for the pure Cahn-Hilliard equation with dynamic boundary conditions}, Adv. Nonlinear Anal., 4 (2015), no. 4, 311-325. \bibitem{Russ3} J. Douglas, T.F. Russel; \emph{Numerical methods for convection dominated diffusion problems based on combining the method of characteristics with finite element methods or finite difference method}, SIAM, J. Numer. Anal., 19 (1982), 871-885. \bibitem{113} P. Ghadimi, M. Y. Fard, A. Dashtimanesh; \emph{Application of an Iterative High Order Difference Scheme Along With an Explicit System Solver for Solution of Stream Function-Vorticity Form of Navier-Stokes Equations}, Journal of fluids engineering-transactions of the ASME, 135(4), 2013. \bibitem{13} V. Girault, P. A. Raviart; \emph{Finite element methods for Navier-Stokes equations}, Theory and Algorithms, (1986) Springer Verlag, Berlin. \bibitem{14} V. Girault, J. Giroire, A. Sequeira; \emph{A stream function-vorticity variational formulation for the exterior Stokes problem in weighted Sobolev spaces}, Math. Meth. in the Applied Sciences, 15(5) (1992), 345-363. \bibitem{111} P. Minev, P. N. Vabishchevich; \emph{An operator-splitting scheme for the stream function-vorticity formulation of the unsteady Navier-Stokes equations}, Journal of computational and applied mathematics, 293 (2016), 147-163. \bibitem{17} O. Pironneau; \emph{On the transport diffusion algorithm and its applications to the Navier-Stokes equations}, Numer Math., 38 (1982), 309-332. \bibitem{saa} Y. Saad; \emph{Iterative methods for sparse linear systems}, PWS series in computer science, (1996). \bibitem{333} J. Jaro\v{s}, K. Taka\^si; \emph{Strongly increasing solutions of cyclic systems of second order differential equations with power-type nonlinearities}, Opuscula Math. 35 (2015), no. 1, 47-69. \bibitem{LZL} H. Li, Z. Zhao, Z. Luo; \emph{A space-time continuous finite element method for 2D viscoelastic wave equation}, Bound. Value Probl. 2016, 2016:53, 1-17. \bibitem{HPmpi} HP MPI;; \emph{User's Guide}, Fourth edition, Helwet Packard, Feb 1999. \bibitem{112} W. W. Ren, J. Wu, C. Shu; \emph{A stream function-vorticity formulation-based immersed boundary method and its applications}, International journal for numerical methods in fluids, 70(5) (2012), 627-645. \bibitem{222} Y. Takei; \emph{On the multisummability of WKB solutions of certain singularly perturbed linear ordinary differential equations}, Opuscula, Math. 35 (2015), no. 5, 775-802. \end{thebibliography} \end{document}