The Poisson problem using Nitsche method on two subdomains

The Poisson problem using Nitsche method on two subdomains#

We consider a domain \(\Omega = \Omega_1 \bigcup \Omega_2 = (0,1)^2\), where \(\Omega_1 = (0,\frac{1}{2}) \times (0,1)\) and \(\Omega_2 = (\frac{1}{2}, 1) \times (0,1)\)

We revisit the Poisson equation using the Nitsche method for two subdomains.

\[ \begin{align} - \nabla^2 u = f \quad \text{in $\Omega$}, \quad \quad u = 0 \quad \text{on $\partial \Omega$}, \quad \quad [\![ u ]\!] = 0 \quad \text{on $\mathcal{I}$}, \quad \quad [\![\partial_n u ]\!] = 0 \quad \text{on $\mathcal{I}$}. \end{align} \]

Variational Formulation#

The variational formulation 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} \nabla u \cdot \nabla v ~ d\Omega + \int_{\mathcal{I}} \left( \kappa [\![ u ]\!] ~ [\![ v ]\!] - [\![ u ]\!] ~ \partial_n v - \partial_n u ~ [\![ v ]\!] \right) ~ d\mathcal{I} \),

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

Formal Model#

from sympde.expr     import BilinearForm, LinearForm, integral, Norm
from sympde.expr     import find, EssentialBC
from sympde.topology import ScalarFunctionSpace, Square, Domain, element_of
from sympde.calculus import grad, dot
from sympde.calculus import jump, avg, minus, plus, Dn
from sympde.core     import Constant
from sympde.topology import NormalVector

from sympy import pi, sin

from psydac.api.discretization import discretize

# ... create a domain as the union of two subdomains
A = Square('A',bounds1=(0, 0.5), bounds2=(0, 1))
B = Square('B',bounds1=(0.5, 1.), bounds2=(0, 1))

connectivity = [((0,0,1),(1,0,-1))]
subdomains = [A,B]
domain = Domain.join(subdomains, connectivity, 'domain')
# ...

# one sided approximation of the normal flux on the interface
nn   = NormalVector('nn')
#Dn_I = lambda u: 0.5*(dot(grad(minus(u)), nn) + dot(grad(plus(u)), nn))
Dn_I = lambda u: 0.5*(plus(Dn(u)) + minus(Dn(u)))

# internal interafaces of the domain
I = domain.interfaces

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

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(u),grad(v)))            
                 + integral(I, kappa * jump(u)*jump(v) - Dn_I(u)*jump(v) - jump(u)*Dn_I(v))
                )

# linear form
f = 2*pi**2*sin(pi*x)*sin(pi*y)
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(kappa=1e3)

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)

0.0002193529744905818

Remarks#

There should be a problem in sympde, since computing the TerminalExpr of the bilinear form gives

>>> from sympde.expr import TerminalExpr
>>> expr = TerminalExpr(a, domain)
>>> print(expr)

(InterfaceExpression(A|B, 
                     -kappa*MinusInterfaceOperator(u)*PlusInterfaceOperator(v) 
                     + (0.5*dx1(MinusInterfaceOperator(u))*nn[0] + 0.5*dx2(MinusInterfaceOperator(u))*nn[1])*PlusInterfaceOperator(v) 
                     + (-0.5*dx1(PlusInterfaceOperator(v))*nn[0] - 0.5*dx2(PlusInterfaceOperator(v))*nn[1])*MinusInterfaceOperator(u)
                    ), 
 InterfaceExpression(A|B, 
                     -kappa*MinusInterfaceOperator(v)*PlusInterfaceOperator(u) 
                     + (0.5*dx1(MinusInterfaceOperator(v))*nn[0] + 0.5*dx2(MinusInterfaceOperator(v))*nn[1])*PlusInterfaceOperator(u) 
                     + (-0.5*dx1(PlusInterfaceOperator(u))*nn[0] - 0.5*dx2(PlusInterfaceOperator(u))*nn[1])*MinusInterfaceOperator(v)
                    ), 
 DomainExpression(A, 
                  dx1(u)*dx1(v) + dx2(u)*dx2(v)
                 ), 
 DomainExpression(B, 
                  dx1(u)*dx1(v) + dx2(u)*dx2(v)
                 ), 
 BoundaryExpression(A_\Gamma_2, 
                    kappa*u*v 
                    + u*(-0.5*dx1(v)*nn[0] - 0.5*dx2(v)*nn[1]) 
                    + v*(-0.5*dx1(u)*nn[0] - 0.5*dx2(u)*nn[1])
                   ), 
 BoundaryExpression(B_\Gamma_1, 
                    kappa*u*v 
                    + u*(-0.5*dx1(v)*nn[0] - 0.5*dx2(v)*nn[1]) 
                    + v*(-0.5*dx1(u)*nn[0] - 0.5*dx2(u)*nn[1])
                   )
)

The BoundaryExpression terms should not appear in the TerminalExpr.

Contents