The Poisson equation with weak imposition of Dirichlet conditions#
In this example, we consider the following Poisson problem
\[
\begin{align}
- \nabla^2 u = f \quad \text{in $\Omega$}, \quad \quad
u = g \quad \text{on $\partial \Omega$}.
\end{align}
\]
The variational formulation#
\[
\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} \nabla u \cdot \nabla v ~ d\Omega + \int_{\partial \Omega} \left( \kappa u v - u \partial_n v - v \partial_n u \right)~d\partial\Omega \),
\(l(v) := \int_{\Omega} f v ~ d\Omega + \int_{\partial\Omega} \left( \kappa g v - g \partial_n v \right) ~ d\Gamma\).
Formal Model#
from sympde.expr import BilinearForm, LinearForm, integral, Norm
from sympde.expr import find, EssentialBC
from sympde.topology import ScalarFunctionSpace, Square, element_of
from sympde.calculus import grad, dot, Dn
from sympde.core import Constant
from psydac.api.discretization import discretize
from sympy import pi, sin
kappa = Constant('kappa', is_real=True)
domain = Square()
V = ScalarFunctionSpace('V', domain)
x,y = domain.coordinates
u,v = [element_of(V, name=i) for i in ['u', 'v']]
# bilinear form
a = BilinearForm((u,v), integral(domain, dot(grad(v), grad(u))))
an = BilinearForm((u,v), integral(domain.boundary, kappa*u*v - u*Dn(v) - v*Dn(u)))
# linear form
ue = sin(pi*x)*sin(pi*y)
g = ue
f = 2*pi**2*sin(pi*x)*sin(pi*y)
l = LinearForm(v, integral(domain, f*v))
ln = LinearForm(v, integral(domain.boundary, kappa*g*v - g*Dn(v)))
# Variational problem
equation = find(u, forall=v, lhs=a(u,v) + an(u,v), rhs=l(v) + ln(v))
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])
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[3], line 8
5 Vh = discretize(V, domain_h, degree=degree)
7 # Discretize equation
----> 8 equation_h = discretize(equation, domain_h, [Vh, Vh])
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/discretization.py:593, in discretize(a, *args, **kwargs)
590 return DiscreteFunctional(a, kernel_expr, *args, **kwargs)
592 elif isinstance(a, sym_Equation):
--> 593 return DiscreteEquation(a, *args, **kwargs)
595 elif isinstance(a, BasicFunctionSpace):
596 return discretize_space(a, *args, **kwargs)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/equation.py:106, in DiscreteEquation.__init__(self, expr, *args, **kwargs)
103 test_space = trial_test[1]
104 # ...
--> 106 self._lhs = discretize(expr.lhs, domain, trial_test, **kwargs)
107 # ...
109 self._rhs = discretize(expr.rhs, domain, test_space, **kwargs)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/discretization.py:577, in discretize(a, *args, **kwargs)
574 kernel_expr = TerminalExpr(a, domain)
576 if len(kernel_expr) > 1:
--> 577 return DiscreteSumForm(a, kernel_expr, *args, **kwargs)
579 # TODO uncomment when the SesquilinearForm subclass of bilinearForm is create in SymPDE
580 # if isinstance(a, sym_SesquilinearForm):
581 # return DiscreteSesquilinearForm(a, kernel_expr, *args, **kwargs)
583 if isinstance(a, sym_BilinearForm):
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/fem.py:1701, in DiscreteSumForm.__init__(self, a, kernel_expr, *args, **kwargs)
1699 elif isinstance(a, sym_BilinearForm):
1700 kwargs['update_ghost_regions'] = False
-> 1701 ah = DiscreteBilinearForm(a, e, *args, assembly_backend=backend, **kwargs)
1702 kwargs['matrix'] = ah._matrix
1703 operator = ah._matrix
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/fem.py:361, in DiscreteBilinearForm.__init__(self, expr, kernel_expr, domain_h, spaces, nquads, matrix, update_ghost_regions, backend, linalg_backend, assembly_backend, symbolic_mapping)
357 linalg_backend = backend or linalg_backend
359 # BasicDiscrete generates the assembly code and sets the following attributes that are used afterwards:
360 # self._func, self._free_args, self._max_nderiv and self._backend
--> 361 BasicDiscrete.__init__(self, expr, kernel_expr, comm=comm, root=0, discrete_space=discrete_space,
362 nquads=nquads, is_rational_mapping=is_rational_mapping, mapping=symbolic_mapping,
363 mapping_space=mapping_space, num_threads=self._num_threads, backend=assembly_backend)
365 #... Handle the special case where the current MPI process does not need to do anything
366 if isinstance(target, (Boundary, Interface)):
367
368 # If process does not own the boundary or interface, do not assemble anything
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/basic.py:298, in BasicDiscrete.__init__(self, expr, kernel_expr, folder, comm, root, discrete_space, nquads, is_rational_mapping, mapping, mapping_space, num_threads, backend)
294 def __init__(self, expr, kernel_expr, *, folder=None, comm=None, root=None, discrete_space=None,
295 nquads=None, is_rational_mapping=None, mapping=None,
296 mapping_space=None, num_threads=None, backend=None):
--> 298 BasicCodeGen.__init__(self, expr, folder=folder, comm=comm, root=root, discrete_space=discrete_space,
299 kernel_expr=kernel_expr, nquads=nquads, is_rational_mapping=is_rational_mapping,
300 mapping=mapping, mapping_space=mapping_space, num_threads=num_threads, backend=backend)
301 # ...
302 self._kernel_expr = kernel_expr
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/basic.py:146, in BasicCodeGen.__init__(self, expr, folder, comm, root, discrete_space, kernel_expr, nquads, is_rational_mapping, mapping, mapping_space, num_threads, backend)
134 # ...
135
136 # ... when using user defined functions, there must be passed as
(...) 142 # # TODO raise appropriate error message
143 # raise ValueError('can not find {} implementation'.format(f))
145 if ast:
--> 146 self._save_code(self._generate_code(), backend=self.backend['name'])
148 if comm is not None and comm.size>1: comm.Barrier()
149 # compile code
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/basic.py:242, in BasicCodeGen._generate_code(self)
233 psydac_ast = self.ast
235 parser_settings = {
236 'dim' : psydac_ast.dim,
237 'nderiv' : psydac_ast.nderiv,
238 'mapping': psydac_ast.mapping,
239 'target' : psydac_ast.domain
240 }
--> 242 pyccel_ast = parse(psydac_ast.expr, settings=parser_settings, backend=self.backend)
243 python_code = pycode(pyccel_ast)
245 return python_code
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:128, in parse(expr, settings, backend)
108 """
109 This function takes a Psydac Ast and returns a Pyccel Ast
110
(...) 125
126 """
127 psy_parser = Parser(settings, backend)
--> 128 ast = psy_parser.doit(expr)
129 return ast
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:197, in Parser.doit(self, expr, **settings)
196 def doit(self, expr, **settings):
--> 197 return self._visit(expr, **settings)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:528, in Parser._visit_DefNode(self, expr, **kwargs)
525 if thread_args:
526 arguments += flatten([self._visit(i, **kwargs) for i in thread_args])
--> 528 body = flatten(tuple(self._visit(i, **kwargs) for i in expr.body))
530 inits = []
531 for k,i in self.shapes.items():
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:528, in <genexpr>(.0)
525 if thread_args:
526 arguments += flatten([self._visit(i, **kwargs) for i in thread_args])
--> 528 body = flatten(tuple(self._visit(i, **kwargs) for i in expr.body))
530 inits = []
531 for k,i in self.shapes.items():
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1071, in Parser._visit_Reduce(self, expr, **kwargs)
1067 stmts = list(loop.stmts) + [Reduction(op, rhs, lhs)]
1068 loop = Loop(loop.iterable, loop.index, stmts=stmts, mask=loop.mask, parallel=parallel,
1069 default=default, shared=shared, private=private,
1070 firstprivate=firstprivate, lastprivate=lastprivate, reduction=reduction)
-> 1071 return self._visit(loop, **kwargs)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1992, in Parser._visit_Loop(self, expr, **kwargs)
1986 geo_stmts = self._visit(geo_stmts, **kwargs)
1987 # ...
1988
1989 # ...
1990 # visit loop statements
-> 1992 stmts = self._visit(expr.stmts, **kwargs)
1993 stmts = flatten(stmts)
1995 # update with product statements if available
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:394, in Parser._visit_Tuple(self, expr, **kwargs)
393 def _visit_Tuple(self, expr, **kwargs):
--> 394 args = [self._visit(i) for i in expr]
395 return Tuple(*args)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:404, in Parser._visit_Block(self, expr, **kwargs)
403 def _visit_Block(self, expr, **kwargs):
--> 404 body = [self._visit(i) for i in expr.body]
405 body = flatten(body)
406 if len(body) == 1:
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1992, in Parser._visit_Loop(self, expr, **kwargs)
1986 geo_stmts = self._visit(geo_stmts, **kwargs)
1987 # ...
1988
1989 # ...
1990 # visit loop statements
-> 1992 stmts = self._visit(expr.stmts, **kwargs)
1993 stmts = flatten(stmts)
1995 # update with product statements if available
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:394, in Parser._visit_Tuple(self, expr, **kwargs)
393 def _visit_Tuple(self, expr, **kwargs):
--> 394 args = [self._visit(i) for i in expr]
395 return Tuple(*args)
[... skipping similar frames: Parser._visit at line 249 (1 times)]
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1992, in Parser._visit_Loop(self, expr, **kwargs)
1986 geo_stmts = self._visit(geo_stmts, **kwargs)
1987 # ...
1988
1989 # ...
1990 # visit loop statements
-> 1992 stmts = self._visit(expr.stmts, **kwargs)
1993 stmts = flatten(stmts)
1995 # update with product statements if available
[... skipping similar frames: Parser._visit at line 249 (1 times)]
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:394, in Parser._visit_Tuple(self, expr, **kwargs)
393 def _visit_Tuple(self, expr, **kwargs):
--> 394 args = [self._visit(i) for i in expr]
395 return Tuple(*args)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1071, in Parser._visit_Reduce(self, expr, **kwargs)
1067 stmts = list(loop.stmts) + [Reduction(op, rhs, lhs)]
1068 loop = Loop(loop.iterable, loop.index, stmts=stmts, mask=loop.mask, parallel=parallel,
1069 default=default, shared=shared, private=private,
1070 firstprivate=firstprivate, lastprivate=lastprivate, reduction=reduction)
-> 1071 return self._visit(loop, **kwargs)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1992, in Parser._visit_Loop(self, expr, **kwargs)
1986 geo_stmts = self._visit(geo_stmts, **kwargs)
1987 # ...
1988
1989 # ...
1990 # visit loop statements
-> 1992 stmts = self._visit(expr.stmts, **kwargs)
1993 stmts = flatten(stmts)
1995 # update with product statements if available
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:394, in Parser._visit_Tuple(self, expr, **kwargs)
393 def _visit_Tuple(self, expr, **kwargs):
--> 394 args = [self._visit(i) for i in expr]
395 return Tuple(*args)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1165, in Parser._visit_Reduction(self, expr, **kwargs)
1162 if not( lhs is None ):
1163 lhs = self._visit(lhs)
-> 1165 return self._visit(expr, op=op, lhs=lhs)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:249, in Parser._visit(self, expr, **settings)
247 annotation_method = '_visit_' + cls.__name__
248 if hasattr(self, annotation_method):
--> 249 return getattr(self, annotation_method)(expr, **settings)
250 # Unknown object, we raise an error.
251 raise NotImplementedError('{}'.format(type(expr)))
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/psydac/api/ast/parser.py:1231, in Parser._visit_ComputeKernelExpr(self, expr, op, lhs, **kwargs)
1229 # Create a new name for the temporaries used in each patch
1230 name = get_name(lhs)
-> 1231 temps, rhs = cse_main.cse(rhs, symbols=cse_main.numbered_symbols(prefix=f'temp{name}'))
1233 normal_vec_stmts = []
1234 normal_vectors = expr.expr.atoms(NormalVector)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympy/simplify/cse_main.py:824, in cse(exprs, symbols, optimizations, postprocess, order, ignore, list)
821 symbols = iter(symbols)
823 # Find other optimization opportunities.
--> 824 opt_subs = opt_cse(reduced_exprs, order)
826 # Main CSE algorithm.
827 replacements, reduced_exprs = tree_cse(reduced_exprs, symbols, opt_subs,
828 order, ignore)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympy/simplify/cse_main.py:529, in opt_cse(exprs, order)
527 for e in exprs:
528 if isinstance(e, (Basic, Unevaluated)):
--> 529 _find_opts(e)
531 # split muls into commutative
532 commutative_muls = OrderedSet()
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympy/simplify/cse_main.py:507, in opt_cse.<locals>._find_opts(expr)
504 return expr
505 seen_subexp.add(expr)
--> 507 list(map(_find_opts, expr.args))
509 if _coeff_isneg(expr):
510 neg_expr = -expr
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympy/simplify/cse_main.py:507, in opt_cse.<locals>._find_opts(expr)
504 return expr
505 seen_subexp.add(expr)
--> 507 list(map(_find_opts, expr.args))
509 if _coeff_isneg(expr):
510 neg_expr = -expr
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympy/simplify/cse_main.py:500, in opt_cse.<locals>._find_opts(expr)
497 return
499 if iterable(expr):
--> 500 list(map(_find_opts, expr))
501 return
503 if expr in seen_subexp:
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympde/calculus/core.py:214, in DiffOperator.__getitem__(self, indices, **kw_args)
212 return Indexed(self, *indices, **kw_args)
213 else:
--> 214 return Indexed(self, indices, **kw_args)
File /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympy/tensor/indexed.py:155, in Indexed.__new__(cls, base, *args, **kw_args)
153 if isinstance(base, (str, Symbol)):
154 base = IndexedBase(base)
--> 155 elif not hasattr(base, '__getitem__') and not isinstance(base, IndexedBase):
156 raise TypeError(filldedent("""
157 The base can only be replaced with a string, Symbol,
158 IndexedBase or an object with a method for getting
159 items (i.e. an object with a `__getitem__` method).
160 """))
161 args = list(map(sympify, args))
KeyboardInterrupt:
Solving the PDE#
equation_h.set_solver('gmres', info=False, tol=1e-8)
uh = equation_h.solve(kappa=1e3)
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)