1
%----------------------------------------------------
2
%\section{Non-homogeneous Neumann boundary conditions for the Laplace operator}
3
%----------------------------------------------------
4
\label{sec-neumann-laplace}
8
In this chapter we study how to solve a ill-posed problem
9
with a solution defined up to a constant.
11
\subsection*{Formulation}
12
Let $f \in L^{2}(\Omega)$ and
13
$g \in H^{\frac{1}{2}}(\partial \Omega)$ satisfying
14
the following compatibility condition:
16
\int_\Omega f \, {\rm d}x
18
\int_{\partial\Omega} g \, {\rm d}s
23
$(P_5)_h$: {\it find $u$, defined in $\Omega$ such that:}
25
-\Delta u &=& f \ {\rm in}\ \Omega \\
26
\Frac{\partial u}{ \partial n} &=& g \ {\rm on}\ \partial \Omega
28
\cindex{algorithm!conjugate gradient}
29
\cindex{algorithm!minres}
30
Since this problem only involves the derivatives of $u$,
31
it is clear that its solution is never unique~\cite[p.~11]{girault-raviart-86}.
32
A discrete version of this problem could be solved iteratively
33
by the conjugate gradient or the MINRES algorithm~\cite{PaiSau-1975}.
34
In order to solve it by a direct method,
35
we turn the difficulty by seeking $u$ in the following space
37
V = \{ v\in H^1(\Omega);\ \ b(v,1) = 0 \}
41
b(v,\mu) = \int_\Omega v \, {\rm d}x,
42
\ \ \forall v\in L^2(\Omega),
43
\forall \mu\in \mathbb{R}
45
The variational formulation of this problem writes:
47
$(VF_5)$: {\it find $u \in V$ such that:}
49
a(u,v) = l(v), \ \forall v \in V
53
a(u,v) &=& \int_\Omega \nabla u . \nabla v \, dx \\
54
l(v) &=& m(f,v) + m_b(g, v) \\
55
m(f,v) &=& \int_\Omega f v \, dx \\
56
m_b(g,v) &=& \int_{\partial \Omega} g v \, ds
58
Since the direct discretization of the space $V$ is not
60
the constraint $b(u,1) = 0$ is enforced by a Lagrange multiplier
61
\cindex{Lagrange!multiplier}%
62
$\lambda\in \mathbb{R}$.
63
Let us introduce the Lagrangian, defined
64
for all $v\in H^1(\Omega)$ and $\mu\in\mathbb{R}$ by:
66
L(v,\mu) = \Frac{1}{2} a(v,v) + b(v,\mu) - l(v)
68
The saddle point $(u,\lambda)\in H^1(\Omega)\times\mathbb{R}$
69
of this Lagrangian is characterized as the unique solution
73
&=& l(v), \ \ \forall v\in H^ 1(\Omega) \\
74
b(u,\mu) \phantom{+ b(v,\lambda)}
75
&=& 0, \ \ \ \ \ \forall \mu\in\mathbb{R}
77
It is clear that if $(u,\lambda)$ is solution of this problem,
78
then $u\in V$ and $u$ is a solution of $(VF_5)$.
79
Conversely, let $u\in V$ the solution of $(VF_5)$.
80
Choosing $v=v_0$ where $v_0(x)=1$, $\forall x\in\Omega$
83
\lambda\,{\rm meas}(\Omega)=l(v_0)
85
From the definition of $l(.)$ and the compatibility
86
condition between the data $f$ and $g$, we get $\lambda=0$.
87
Notice that the saddle point problem extends to the case
88
when $f$ and $g$ does not satisfies the compatibility condition,
89
and in that case $\lambda=l(v_0)/{\rm meas}(\Omega)$.
91
\subsection*{Approximation}
92
\cindex{Lagrange!interpolation}
93
As usual, we introduce a mesh ${\cal T}_h$ of $\Omega$
94
and the finite dimensional space $X_h$:
96
X_h = \{ v \in H^1(\Omega); \
98
\forall K \in {\cal T}_h \}
100
The approximate problem writes:\\
101
{\it $(VF_5)_h$: find $(u_h,\lambda_h)\in X_h\times \mathbb{R}$ such that:}
103
a(u_h,v) + b(v,\lambda_h)
104
&=& l_h(v), \ \ \forall v\in X_h \\
105
b(u_h,\mu) \phantom{+ b(v,\lambda)}
106
&=& 0, \ \ \ \ \ \forall \mu\in\mathbb{R}
110
l_h(v) = m(\Pi_h f,v_h) + m_b(\pi_h g, v_h)
113
\subsection*{File \file{neumann-laplace.cc}}
114
\exindex{neumann-laplace.cc}
115
\myexample{neumann-laplace.cc}
117
\subsection*{Comments}
118
Let $\Omega \subset R^d$, $d=1,2,3$.
119
We choose $f(x) = 1$ and $g(x) = -1/(2d)$.
120
This example is convenient, since the exact solution is known:
122
u(x) = - \Frac{1}{12} + \Frac{1}{2d} \sum_{i=1}^d x_i(1-x_i)
124
The code looks like the previous ones.
125
Let us comment the changes.
126
The discrete bilinear form $b$ is computed as $b_h\in X_h$
127
that interprets as a linear application from $X_h$ to $\mathbb{R}$:
128
$b_h(v_h)=m(v_h,1)$. Thus $b_h$ is computed as
129
\begin{lstlisting}[numbers=none,frame=none]
130
field b = m*field(Xh,1.0);
132
where the discrete bilinear form $m$ is identified to its matrix
133
and {\tt field(Xh,1.0)} is the constant vector equal to $1$.
138
\left( \begin{array}{cc}
139
{\tt a.uu} & {\tt trans(b.u)} \\
145
\left( \begin{array}{c}
152
\left( \begin{array}{c}
157
The problem admits the following matrix form:
159
\mathcal{A} \ \mathcal{U} = \mathcal{B}
163
\cindex{matrix!concatenation}
164
The matrix and right-hand side of the system are assembled by concatenation:
165
\begin{lstlisting}[numbers=none,frame=none]
166
csr<Float> A = {{ a.uu, b.u},
168
vec<Float> B = { lh.u, 0 };
171
{\color{rheoleftypcolor} \code{csr}}
173
{\color{rheoleftypcolor} \code{vec}}
174
are respectively the matrix and vector classes.
175
The {\color{rheoleftypcolor} \code{csr}} is the abbreviation
176
of {\em compressed sparse row}, a sparse matrix compression standard format.
178
Notice that the matrix $\mathcal{A}$ is symmetric and non-singular, but indefinite~:
179
it admits eigenvalues that are either strictly positive or strictly negative.
180
Then, the \code{uh.u} vector is extracted from the \code{U} one:
181
\begin{lstlisting}[numbers=none,frame=none]
182
uh.u = U [range(0,uh.u.size())];
184
The extraction of \code{lambda} from \code{U} is more technical in a distributed
185
environment. In a sequential one, since it is stored after the \code{uh.u} values,
186
it could be simply written as:
187
\begin{lstlisting}[numbers=none,frame=none]
188
Float lambda = U [uh.u.size()];
190
\cindex{distributed computation}%
191
\cindex{parallel computation}%
193
\toindex{library!boost}%
194
\clindex{communicator}%
195
In a distributed environment, \code{lambda} is stored in \code{U} on
196
the last processor, identified by \code{U.comm().size()-1}.
197
Here \code{U.comm()} denotes the \code{communicator}, from the \code{boost::mpi}
198
library and \code{U.comm().size()} is the number of processors in use,
199
e.g. as specified by the \code{mpirun} command.
200
On this last processor, the array \code{U} has size equal to
201
\code{uh.u.size()+1} and \code{lambda} is stored in \code{U[uh.u.size()]}.
202
On the others processors, the array \code{U} has size equal to
203
\code{uh.u.size()} and lambda is not available.
204
The following statement extract \code{lambda} on the last processor
205
and set it to zero on the others:
206
\begin{lstlisting}[numbers=none,frame=none]
207
Float lambda = (U.size() == uh.u.size()+1) ? U [uh.u.size()] : 0;
209
Then, the value of \code{lambda} is broadcasted on the others processors:
210
\begin{lstlisting}[numbers=none,frame=none]
211
mpi::broadcast (U.comm(), lambda, U.comm().size() - 1);
213
The preprocessing guards \code{\#idef}$\ldots$\code{\#endif} assure that
214
this example compile also when the library is not installed with the MPI support.
216
Finally, the statement
217
\begin{lstlisting}[numbers=none,frame=none]
218
dout << catchmark("u") << uh
219
<< catchmark("lambda") << lambda << endl;
221
writes the solution $(u_h,\lambda)$.
222
The \code{catchmark} function writes marks together with the solution in the output stream.
223
These marks are suitable for a future reading with the same format, as:
224
\begin{lstlisting}[numbers=none,frame=none]
225
din >> catchmark("u") >> uh
226
>> catchmark("lambda") >> lambda;
228
This is usefull for post-treatment, visualization and error analysis.
230
\subsection{How to run the program}
235
mkgeo_grid -t 10 > square.geo
236
./neumann-laplace square P1 | field -
added
removed
doc/pusrman/neumann-laplace.tex
