Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Domain Decomposition Methods
Author
By Abdoulaye Samake
Vincent Chabannes
Christophe Prud’homme




Introduction

In mathematics, numerical analysis, and numerical partial differential equations, domain decomposition methods solve a boundary value problem by splitting it into smaller boundary value problems on subdomains and iterating to coordinate the solution between the adjacent subdomains. A corse problem with one or fiew unknows per subdomain is used to further coordinate the solution between the subdomains globally.

A 1D Model

We consider the following laplacian boundary value problem

\begin{equation} \left \{ \begin{aligned} & -u"(x) = f(x) \quad \text{in} \quad ]0,1[ \\ & u(0) =\alpha, ~ u(1) = \beta \end{aligned} \right. \label{eq:30} \end{equation}

where $\alpha, \beta \in \mathbb R.$

Schwartz Algorithms

The schwartz overlapping multiplicative algorithm with dirichlet interface conditions for this problem at $n^{th}$ iteration is given by

\begin{equation} \label{eq:31} \left \{ \begin{aligned} -u_1"^n(x) & = f(x) \quad \text{in} \quad ]0,b[ \\ u_1^n(0) & = \alpha \\ u_1^n(b) & = u_2^{n-1}(b) \end{aligned} \right. \qquad \text{and} \qquad \left \{ \begin{aligned} -u_2"^n(x) & = f(x) \quad \text{in} \quad ]a,1[ \\ u_2^n(1) & = \beta \\ u_2^n(a) & = u_1^n(a) \end{aligned} \right. \end{equation}

where $ n \in \mathbb N^*, a, b \in \mathbb R $ and $a < b$.
Let $e_i^n = u_i^n-u~(i=1,2)$, the error at $n^{th}$ iteration relative to the exact solution, the convergence rate is given by

\begin{equation} \rho = \frac{\vert e_1^n \vert}{\vert e_1^{n-1} \vert} = \frac{a}{b}\frac{1-b}{1-a} = \frac{\vert e_2^n \vert}{\vert e_2^{n-1} \vert} . \label{eq:32} \end{equation}

Variational formulations

find $u$ such that

\begin{equation*} \int_0^b u_1'v' = \int_0^b fv \quad \forall v \qquad \text{in the first subdomain} ~\Omega_1 = ]0,b[ \end{equation*}

\begin{equation*} \int_a^1 u_2'v' = \int_a^1 fv \quad \forall v \qquad \text{in the second subdomain} ~ \Omega_2 = ]a,1[ \end{equation*}

A 2 domain overlapping Schwartz method in 2D and 3D

We consider the following laplacian boundary value problem

\begin{equation} \left \{ \begin{aligned} -\Delta u & = f \quad \text{in} \quad \Omega \\ u & = g \quad \text{on} \quad \partial\Omega \end{aligned} \right. \label{eq:33} \end{equation}

where $\Omega \subset \mathbb R^d, d=2,3$ and $g$ is the dirichlet boundary value.

Schwartz algorithms

The schwartz overlapping multiplicative algorithm with dirichlet interface conditions for this problem on two subdomains $\Omega_1$ and $\Omega_2$ at $n^{th}$ iteration is given by

\begin{equation} \label{eq:34} \left \{ \begin{aligned} -\Delta u_1^n & = f \quad \qquad \text{in} \quad \Omega_1 \\ u_1^n & = g \quad \qquad \text{on} \quad \partial \Omega_1^{ext}\\ u_1^n & = u_2^{n-1} \quad ~~ \text{on} \quad \Gamma_1 \end{aligned} \right. \qquad \text{and} \qquad \left \{ \begin{aligned} - \Delta u_2^n & = f \quad \qquad \text{in} \quad \Omega_2 \\ u_2^n & = g \quad \qquad \text{on} \quad \partial \Omega_2^{ext}\\ u_2^n & = u_1^n \qquad~~ \text{on} \quad \Gamma_2 \end{aligned} \right. \end{equation}

Variational formulations

\begin{equation*} \begin{aligned} \int_{\Omega_i} \nabla u_i \cdot \nabla v = \int_{\Omega_i} fv \quad \forall~ v,~i=1,2. \end{aligned} \end{equation*}

Implementation

// Implementation of the local problem
template<Expr>
void
localProblem(element_type& u, Expr expr)
{
// Assembly of the right hand side $ \mathlarger \int_\Omega fv $
auto F = M_backend->newVector(Xh);
form1( _test=Xh,_vector=F, _init=true ) =
integrate( elements(mesh), f*id(v) );
F->close();
// Assembly of the left hand side $ \mathlarger \int_\Omega \nabla u \cdot \nabla v$
auto A = M_backend->newMatrix( Xh, Xh );
form2( _test=Xh, _trial=Xh, _matrix=A, _init=true ) =
integrate( elements(mesh), gradt(u)*trans(grad(v)) );
A->close();
// Apply the dirichlet boundary conditions
form2( Xh, Xh, A ) +=
on( markedfaces(mesh, "Dirichlet") ,u,F,g);
// Apply the dirichlet interface conditions
form2( Xh, Xh, A ) +=
on( markedfaces(mesh, "Interface") ,u,F,expr);
// solve the linear system $ A u = F $
M_backend->solve(_matrix=A, _solution=u, _rhs=F );
}
unsigned int cpt = 0;
double tolerance = $1e-8$;
double maxIterations = 20;
double l2erroru$_1$ = 1.;
double l2erroru$_2$ = 1;
// Iteration loop
while( (l2erroru$_1$ +l2erroru$_2$) > tolerance && cpt <= maxIterations)
{
// call the localProblem on the first subdomain $\Omega_1$
localProblem(u$_1$, idv(u$_2$));
// call the localProblem on the first subdomain $\Omega_2$
localProblem(u$_2$, idv(u$_1$));
// compute L2 errors on each subdomain
L2erroru$_1$ = l2Error(u$_1$);
L2erroru$_2$ = l2Error(u$_2$);
// increment the cunter
++cpt;
}

Numerical results in 2D case

The numerical results presented in the following table correspond to the partition of the global domain $\Omega$ in two subdomains $\Omega_1$ and $\Omega_2$ and the following configuration:

  • $ g(x,y) = \sin(\pi x)\cos(\pi y)$ : the exact solution
  • $f(x,y) = 2\pi^2g$ : the right hand side of the equation
  • $\mathbb P_2$ approximation : the lagrange polynomial order
  • hsize $= 0.02$ : the mesh size
  • tol $=1e-9$ : the tolerance
dd2dgeometry.png
dd2dmesh.png
Two overlapping subdomains
Two overlapping meshes
Geometry
Nomber of iterations$\mathbf {\| u_1-u_{ex}\|_{L_2} }$$\mathbf{\| u_2-u_{ex}\|_{L_2}}$
11 2.52e-8 2.16e-8

Numerical Solutions in 2D Case

iter_1.png
iter_10.png
First Iteration
$ 10^{th}$ Iteration
Isovalues of Solution in 2D

Computing the eigenmodes of the Dirichlet to Neumann operator

Problem description and variational formulation

We consider at the continuous level the Dirichlet-to-Neumann(DtN) map on $\Omega$, denoted by DtN $_{\Omega}$. Let $u: \Gamma \longmapsto \mathbb R, $

\begin{equation*} \label{eq:35} \text{DtN}_{\Omega}(u) = \kappa \frac{\partial v}{ n} \Big |_{\Gamma} \end{equation*} where $v$ satisfies \begin{equation} \label{eq:36} \left\{ \begin{aligned} & \mathcal L(v):= (\eta - \text{div}(\kappa \nabla))v = 0 & \text{dans} \quad \Omega,\\ & v = u & \text{sur} \quad \Gamma \end{aligned} \right. \end{equation*}

where $\Omega$ is a bounded domain of $\mathbb R^d$ (d=2 or 3), and $\Gamma$ it border, $\kappa$ is a positive diffusion function which can be discontinuous, and $\eta \geq 0$. The eigenmodes of the Dirichlet-to-Neumann operator are solutions of the following eigenvalues problem

\begin{equation} \label{eq:37} \text{DtN}_{\Omega}(u) = \lambda \kappa u \end{equation}

To obtain the discrete form of the DtN map, we consider the variational form of $(1)$. let's define the bilinear form $a : H^1(\Omega) \times H^1(\Omega) \longrightarrow \mathbb R $,

\begin{equation*} \label{eq:41} a(w,v) := \int_\Omega \eta w v + \kappa \nabla w \cdot \nabla v . \end{equation*}

With a finite element basis $\{ \phi_k \}$, the coefficient matrix of a Neumann boundary value problem in $\Omega$ is

\begin{equation*} \label{eq:42} A_{kl} := \int_\Omega \eta \phi_k \phi_l + \kappa \nabla \phi_k \cdot \nabla \phi_l . \end{equation*}

A variational formulation of the flux reads

\begin{equation*} \label{eq:43} \int_\Gamma \kappa \dfrac{\partial v}{\partial n} \phi_k = \int_\Omega \eta v \phi_k + \kappa \nabla v \cdot \nabla \phi_k \quad \forall~ \phi_k. \end{equation*}

So the variational formulation of the eigenvalue problem $(2)$ reads

\begin{equation} \label{eq:40} \int_\Omega \eta v \phi_k + \kappa \nabla v \cdot \nabla \phi_k = \lambda \int_\Gamma \kappa v \phi_k \quad \forall~ \phi_k. \end{equation}

Let $B$ be the weighted mass matrix

\begin{equation*} \label{eq:44} (B)_{kl} = \int_\Gamma \kappa \phi_k \phi_l \end{equation*}

The compact form of $(3)$ is

\begin{equation} \label{eq:45} Av = \lambda B v \end{equation}

Implementation

Assembly of the right hand side $ B = \mathlarger \int_\Gamma \kappa v w $

auto B = M_backend->newMatrix( Xh, Xh ) ;
form2( _test=Xh, _trial=Xh, _matrix=B, _init=true );
BOOST_FOREACH( int marker, flags )
{
form2( Xh, Xh, B ) +=
integrate( markedfaces(mesh,marker), kappa*idt(u)*id(v) );
}
B->close();

Assembly of the left hand side $ A = \mathlarger \int_\Omega \eta v w + \kappa \nabla v \cdot \nabla w $

auto A = M_backend->newMatrix( Xh, Xh ) ;
form2( _test=Xh, _trial=Xh, _matrix=A, _init=true ) =
integrate( elements(mesh), kappa*gradt(u)*trans(grad(v)) + nu*idt(u)*id(v) );
A->close();
// eigenvalue solver options
int nev = this->vm()["solvereigen-nev"].template as<int>();
int ncv = this->vm()["solvereigen-ncv"].template as<int>();;
// definition of the eigenmodes
SolverEigen<double>::eigenmodes_type modes;

solve the eigenvalue problem $ Av = \lambda B v $

modes=
eigs( _matrixA=A,
_matrixB=B,
_nev=nev,
_ncv=ncv,
_transform=SINVERT,
_spectrum=SMALLEST_MAGNITUDE,
_verbose = true );
}

Numerical solutions

mode-0.png
mode-1.png
mode-2.png
first mode
second mode
third mode
Three eigenmodes

These numerical solutions correspond to the following configuration :

  • $\mathbb P_2$ approximation : the lagrange polynomial order
  • hsize $= 0.02$ : the mesh size
  • $\mu = \kappa = 1.$

Generated on Fri Oct 25 2013 14:24:27 for Feel++ by doxygen 1.8.4