Vector Poisson equation

Vector Poisson equation#

In this example we consider the vector Poisson equation with homogeneous Dirichlet boundary conditions:

\[ \begin{align} - \nabla^2 \mathbf{u} = \mathbf{f} \quad \mbox{in} ~ \Omega, \quad \quad \mathbf{u} = 0 \quad \mbox{on} ~ \partial \Omega. \end{align} \]

The Variational Formulation#

The corresponding variational formulation, using \(\mathbf{H}^1\) formulation, i.e. all components are in \(H^1\), reads

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

where

  • \(V \subset \mathbf{H}_0^1(\Omega)\),

  • \(a(\mathbf{u},\mathbf{v}) := \int_{\Omega} \nabla \mathbf{u} : \nabla \mathbf{v} ~ d\Omega\),

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

Formal Model#

from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC, Norm, SemiNorm
from sympde.topology import VectorFunctionSpace, Cube, element_of
from sympde.calculus import grad, inner, dot

from sympy import pi, sin, Tuple, Matrix

domain = Cube()

V = VectorFunctionSpace('V', domain)

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 , inner(grad(v), grad(u))))

# linear form
f1 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)
f2 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)
f3 = 3*pi**2*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()

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[5], line 1
----> 1 uh = equation_h.solve()

File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/equation.py:198, in DiscreteEquation.solve(self, **kwargs)
    196 def solve(self, **kwargs):
--> 198     self.assemble(**kwargs)
    200     # Free arguments of current equation
    201     free_args = set(self.lhs.free_args + self.rhs.free_args)

File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/equation.py:178, in DiscreteEquation.assemble(self, **kwargs)
    176 # Matrix (left-hand side)
    177 if assemble_lhs:
--> 178     A = self.lhs.assemble(reset=True, **kwargs)
    179     if self.bc:
    180         apply_essential_bc(A, *self.bc)

File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/fem.py:532, in DiscreteBilinearForm.assemble(self, reset, **kwargs)
    529 if reset:
    530     reset_arrays(*self.global_matrices)
--> 532 self._func(*args, *self._threads_args)
    533 if self._matrix and self._update_ghost_regions:
    534     self._matrix.exchange_assembly_data()

File ~/work/IGA-Python/IGA-Python/chapter2/__psydac__/dependencies_lr7xfr33.py:97, in assemble_matrix_lr7xfr33(global_test_basis_v_1, global_test_basis_v_2, global_test_basis_v_3, global_trial_basis_u_1, global_trial_basis_u_2, global_trial_basis_u_3, global_span_v_1, global_span_v_2, global_span_v_3, global_x1, global_x2, global_x3, test_v_p1, test_v_p2, test_v_p3, trial_u_p1, trial_u_p2, trial_u_p3, n_element_1, n_element_2, n_element_3, k1, k2, k3, pad1, pad2, pad3, g_mat_u_0_v_0_lr7xfr33, g_mat_u_1_v_1_lr7xfr33, g_mat_u_2_v_2_lr7xfr33)
     95             contribution_v_0_u_0_lr7xfr33 += u_0_x1*v_0_x1 + u_0_x2*v_0_x2 + u_0_x3*v_0_x3
     96             contribution_v_1_u_1_lr7xfr33 += u_1_x1*v_1_x1 + u_1_x2*v_1_x2 + u_1_x3*v_1_x3
---> 97             contribution_v_2_u_2_lr7xfr33 += u_2_x1*v_2_x1 + u_2_x2*v_2_x2 + u_2_x3*v_2_x3
    101 l_mat_u_0_v_0_lr7xfr33[i_basis_1,i_basis_2,i_basis_3,2 - i_basis_1 + j_basis_1,2 - i_basis_2 + j_basis_2,2 - i_basis_3 + j_basis_3] = contribution_v_0_u_0_lr7xfr33
    102 l_mat_u_1_v_1_lr7xfr33[i_basis_1,i_basis_2,i_basis_3,2 - i_basis_1 + j_basis_1,2 - i_basis_2 + j_basis_2,2 - i_basis_3 + j_basis_3] = contribution_v_1_u_1_lr7xfr33

KeyboardInterrupt:

Computing the error norm#

When the analytical solution is available, you might be interested in computing the \(L^2\) norm or \(H^1_0\) semi-norm. SymPDE allows you to do so, by creating the Norm object. In this example, the analytical solution is given by

\[ u_e = \sin(\pi x) \sin(\pi y) \sin(\pi z) \]

Computing the \(L^2\) norm#

ue1 = sin(pi*x)*sin(pi*y)*sin(pi*z)
ue2 = sin(pi*x)*sin(pi*y)*sin(pi*z)
ue3 = sin(pi*x)*sin(pi*y)*sin(pi*z)
ue = Tuple(ue1, ue2, ue3)

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

error = Matrix([u[0]-ue[0], u[1]-ue[1], u[2]-ue[2]])

# create the formal Norm object
l2norm = Norm(error, 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(error, 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(error, 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