Elliptic-div Problem

Elliptic-div Problem#

Let \(\Omega \subset \mathbb{R}^d\) be an open Liptschitz bounded set, and we look for the solution of the following problem

\[\begin{split} \begin{align} \left\{ \begin{array}{rl} - \nabla \nabla \cdot \mathbf{u} + \mu \mathbf{u} &= \mathbf{f}, \quad \Omega \\ \mathbf{u} \times \mathbf{n} &= 0, \quad \partial\Omega \end{array} \right. \end{align} \end{split}\]

where \(\mathbf{f} \in \mathbf{L}^2(\Omega)\), \(\mu \in L^\infty(\Omega)\) and there exists \(\mu_0 > 0\) such that \(\mu \geq \mu_0\) almost everywhere.

The Variational Formulation#

We take the Hilbert space \(V := \mathbf{H}\_0(\mbox{div}, \Omega)\), in which case the variational formulation writes

Find \(\mathbf{u} \in V\) such that

\[ \begin{align} a(\mathbf{u},\mathbf{v}) = l(\mathbf{v}) \quad \forall \mathbf{v} \in V \end{align} \]

where

\[\begin{split} \begin{align} \left\{ \begin{array}{rll} a(\mathbf{u}, \mathbf{v}) &:= \int_{\Omega} \nabla \cdot \mathbf{u} ~ \nabla \cdot \mathbf{v} + \int_{\Omega} \mu \mathbf{u} \cdot \mathbf{v}, & \forall \mathbf{u}, \mathbf{v} \in V \\ l(\mathbf{v}) &:= \int_{\Omega} \mathbf{v} \cdot \mathbf{f}, & \forall \mathbf{v} \in V \end{array} \right. \end{align} \end{split}\]

We recall that in \(\mathbf{H}\_0(\mbox{div}, \Omega)\), the bilinear form \(a\) is equivalent to the inner product and is therefor continuous and coercive. Hence, our abstract theory applies and there exists a unique solution.

Warning

This NB needs to be fixed

Formal Model#

from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC
from sympde.topology import VectorFunctionSpace, Cube, element_of
from sympde.calculus import div, dot
from sympde.core     import Constant

from sympy import pi, sin, Tuple

mu = Constant('mu', is_real=True)

domain = Cube()

V = VectorFunctionSpace('V', domain, kind='Hdiv')

x,y,z = domain.coordinates

u,v = [element_of(V, name=i) for i in ['u', 'v']]

# bilinear form
a = BilinearForm((u,v), integral(domain, div(v)*div(u) + mu * dot(u,v)))

# linear form
f1 = sin(pi*x)*sin(pi*y)*sin(pi*z)
f2 = sin(pi*x)*sin(pi*y)*sin(pi*z)
f3 = sin(pi*x)*sin(pi*y)*sin(pi*z)
f = Tuple(f1, f2, f3)

l = LinearForm(v, integral(domain, dot(f,v)))

# Dirichlet boundary conditions
bc = [EssentialBC(u, 0, domain.boundary)]

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)

Discretization#

We shall need the discretize function from PsyDAC.

from psydac.api.discretization import discretize

degree = [2,2,2]
ncells = [8,8,8]

# Create computational domain from topological domain
domain_h = discretize(domain, ncells=ncells, comm=None)

# Create discrete spline space
Vh = discretize(V, domain_h, degree=degree)

# Discretize equation
equation_h = discretize(equation, domain_h, [Vh, Vh])

Solving the PDE#

#uh = equation_h.solve()

Computing the error norm#

Computing the \(L^2\) norm#

#ue = sin(pi*x)*sin(pi*y)

#u = element_of(V, name='u')

# create the formal Norm object
#l2norm = Norm(u - ue, domain, kind='l2')

# discretize the norm
#l2norm_h = discretize(l2norm, domain_h, Vh)

# assemble the norm
#l2_error = l2norm_h.assemble(u=uh)

# print the result
#print(l2_error)

Computing the \(H^1\) semi-norm#

# create the formal Norm object
#h1norm = SemiNorm(u - ue, domain, kind='h1')

# discretize the norm
#h1norm_h = discretize(h1norm, domain_h, Vh)

# assemble the norm
#h1_error = h1norm_h.assemble(u=uh)

# print the result
#print(h1_error)

Computing the \(H^1\) norm#

# create the formal Norm object
#h1norm = Norm(u - ue, domain, kind='h1')

# discretize the norm
#h1norm_h = discretize(h1norm, domain_h, Vh)

# assemble the norm
#h1_error = h1norm_h.assemble(u=uh)

# print the result
#print(h1_error)
Contents