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)