Elliptic equation in the general form#

We consider here, the following general form of an elliptic partial differential equation,

\[\begin{split} \begin{align} - \nabla \cdot \left( A \nabla u \right) + \mathbf{b} \cdot \nabla u + c u &= f, \quad \Omega \\ u &= 0, \quad \partial \Omega \end{align} \end{split}\]

Variational Formulation#

An \(H^1\)-conforming variational formulation of reads

\[ \begin{align} \text{find $u \in V$ such that} \quad a(u,v) = l(v) \quad \forall v \in V, \end{align} \]

where

  • \(V \subset H^1(\Omega)\),

  • \(a(u,v) := \int_{\Omega} \left( \left( A \nabla u \right) \cdot \nabla v + \left( \mathbf{b} \cdot \nabla u \right) v + cuv \right) ~ d\Omega\),

  • \(l(v) := \int_{\Omega} f v ~ d\Omega\).

Formal Model#

from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC, Norm, SemiNorm
from sympde.topology import ScalarFunctionSpace, Square, element_of
from sympde.calculus import grad, dot, div
from sympde.core import Vector, Matrix

from sympy import pi, sin

from psydac.api.discretization import discretize

domain = Square()

V = ScalarFunctionSpace('V', domain)

x,y = domain.coordinates

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

c = x*y
b = Vector([1e-2, 1e-1], name='b')
A = Matrix([[1,1], [0,1]], name='A')

# bilinear form
expr = dot(grad(v), A * grad(u)) + dot(b, grad(u))*v + c*u*v
a = BilinearForm((u,v), integral(domain, expr))
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 5
      3 from sympde.topology import ScalarFunctionSpace, Square, element_of
      4 from sympde.calculus import grad, dot, div
----> 5 from sympde.core import Vector, Matrix
      7 from sympy import pi, sin
      9 from psydac.api.discretization import discretize

ImportError: cannot import name 'Vector' from 'sympde.core' (/opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympde/core/__init__.py)

Manifactured solution#

# the analytical solution and its rhs
ue = sin(pi * x) * sin(pi * y)

L = lambda u: - div(A*grad(u)) + dot(b,grad(u)) + c*u
f  = L(ue)

Formal Equation#

l = LinearForm(v, integral(domain, 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#

degree = [2,2]
ncells = [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#

equation_h.set_solver('gmres', info=False, tol=1e-8)
uh = equation_h.solve()

Computing the error norm#

Computing the \(L^2\) norm#

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)