Error on SIPG implementation

1

I have been working with Brinkman equantion:

- \mu^{*} \nabla ^{2} \textbf u + \mathbf{u} \frac {\mu} {\textbf{K}} + \nabla p = 0

\nabla \cdot \mathbf{u} = 0 ,

where \mu^{*} and \mu are the Brinkiman viscosity and dynamic viscosity (In this example \mu^{*} and \mu has the same value). \textbf{K} is the permeability tensor, \textbf{u} and p are the velocity and pressure respectively.

I based on my code on this article:
resume or full article

I have implemented some modifications to impose the weak form on pressure boundary. Applying SIPG method we have:

\int_{\Omega}^{}\nabla \mathbf{u} \cdot q dx= 0

\int_{\Omega}^{} \mu^{*} \nabla \textbf u \cdot \nabla v \: dx + \int_{\Omega}^{} \mathbf{u} \frac {\mu} {\textbf{K}} \cdot v \: dx - \int_{\Omega}^{}p \cdot \nabla v \:dx + \int_{\Gamma}^{}n \cdot vp \: ds + Stab= 0

Stab = \mu^{*} \left( \sum_{E \:\epsilon \: \varepsilon }^{} \frac{\alpha }{h}\left \langle \left [ \mathbf{u} \right],\left [ v \right ] \right \rangle_{E} - \left \langle \left \{ \nabla\mathbf{u} \right \},\left [ v \right ] \right \rangle_{E} - \left \langle \left \{ \nabla v \right \},\left [ \mathbf{u} \right ] \right \rangle_{E}\right)

where \left [ \: \right ] and \left \{ \right \} are the jump and average operators, h is the cell size and \alpha is a positive constant on each edge.

MY code is:

from fenics import *


class vugg(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.4, 0.6)) and between(x[0], (0.3, 0.5))


mesh = UnitSquareMesh(100, 100, "crossed")

order = 1
V = FiniteElement("BDM", mesh.ufl_cell(), order)
Q = FiniteElement("DG", mesh.ufl_cell(), order - 1)


Element = V * Q
W = FunctionSpace(mesh, Element)

(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))


left.mark(boundary_parts, 0)
right.mark(boundary_parts, 1)
bottom.mark(boundary_parts, 2)
top.mark(boundary_parts, 3)


obstacle = vugg()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)

ds = Measure("ds", domain=mesh, subdomain_data=boundary_parts)
dx = Measure("dx", domain=mesh, subdomain_data=domains)


mu = Constant(0.001002)
k = Constant(1e-14)
pin = Constant(2)
pout = Constant(0)
alpha = 35
h = CellDiameter(mesh)

f = Constant((0.0, 0.0))
n = FacetNormal(mesh)

stab = (
    mu * (alpha / h) * inner(jump(u), jump(v)) * dS
    - mu * inner(avg(grad(u)), jump(v)) * dS
    - mu * inner(avg(grad(v)), jump(u)) * dS
)

a = (
    mu * inner(grad(u), grad(v)) * dx(1)
    + mu / k * inner(u, v) * dx(0)
    - div(v) * p * dx(1)
    - div(u) * q * dx(0)
    - div(v) * p * dx(0)
    - div(u) * q * dx(1)
    + stab
)


L = (
    inner(f, v) * dx(0)
    + inner(f, v) * dx(1)
    - pin * dot(v, n) * ds(0)
    - pout * dot(v, n) * ds(1)
)

U = Function(W)

solve(a == L, U)

(u1, p1) = U.split()

My code presents the flowing error:

Shapes do not match: <ComponentTensor id=140300829879168> and <Sum id=140300599110272>.
Traceback (most recent call last):
  File "/home/tfk/Desktop/Doutorado/Programas/Brinkman_mono/Brinkman_monophase_BDM_minimun_code.py", line 59, in <module>
    - mu * inner(avg(grad(u)), jump(v)) * dS
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
    return Inner(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
    error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <ComponentTensor id=140300829879168> and <Sum id=140300599110272>.

Any suggestions would be very helpful.

2

ufl's jump operator is defined as seen here. It looks like you’re missing the optional argument to find the jump of the normal flux.

3

I really appreciate your help @nate. But I did not figure out how to deal with this problem, maybe because I am new Fenics user. I have compared my solution to the Biharmonic Equation demo_biharmonic.py and my implementation is similar. I also tried this with no success:

stab = (
    mu * (alpha / h) * inner(jump(u, n), jump(v, n)) * dS
    - mu * inner(avg(grad(u)), jump(v, n)) * dS
    - mu * inner(avg(grad(v)), jump(u, n)) * dS
)

and this:

stab = (
    mu * (alpha / h("+")) * inner(jump(u, n), jump(v, n)) * dS
    - mu * inner(avg(grad(u)), dot(v("+"), n) - dot(v("-"), n)) * dS
    - mu * inner(avg(grad(v)), dot(u("+"), n) - dot(u("-"), n)) * dS
)

Do you have any example or could detail the solution little bit?

Thanks

4

I failed to notice u and v are vector quantities. Sadly ufl doesn’t generalise as far as tensor fluxes. You could use the following which is the ufl representation of

\underline{[\![ \vec{u} ]\!]} = \vec{u}^+ \otimes \vec{n}^+ + \vec{u}^- \otimes \vec{n}^-

import ufl
def tensor_jump(v, n):
    return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")

And replace jump with tensor_jump where appropriate

Also note the quantity h is multi valued on a facet, so consider h = ufl.Min(h("+"), h("-")).

1 Like
5

Thanks @nate , the code is working now

1 Like