Assembling matrix via inner(div(u), div(v)) failed

[thoeschler]

Hi Mikael,

I am trying to solve a linear elasticity problem on a 2D domain. I'm using Legendre polynomials and, thus, I have assembled the weak form using integration by parts, which contains to terms:

1. The inner product of the gradients of the test and trial functions inner(grad(u), grad(v))
2. The product of the divergences of the test and trial functions which I tried to implement using inner(div(u), div(v)).

The first term works fine. The second term inner(div(u), div(v)) raises an error. However, I found a workaround, see my code below. I think it would be convenient to be able to directly use inner(div(u), div(v)).

Best regards, Thilo

from sympy import symbols, cos, sin
from mpi4py import MPI
from shenfun import inner, grad, div, Dx, TestFunction, TrialFunction, Array,\
    Basis, TensorProductSpace, VectorTensorProductSpace, BlockMatrix
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD

family = 'legendre'

# Size of discretization
n = 40
N = [n, n+1]

# Bases and spaces
BX_hom = Basis(N[0], family=family, bc=(0, 0))
BY_hom = Basis(N[1], family=family, bc=(0, 0))

subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T_hom = TensorProductSpace(subcomms, (BX_hom, BY_hom))

V = VectorTensorProductSpace([T_hom, T_hom])

# Function spaces
u = TrialFunction(V)
v = TestFunction(V)

# Compute matrices
# Assemble standard term
A = inner(grad(u), grad(v))

# Assemble (div(u), div(v))-term

# non-working
B = inner(div(u), div(v)) # raises an error

# workaround
C00 = inner(Dx(u[0], 0), Dx(v[0], 0))
C01 = inner(Dx(u[0], 0), Dx(v[1], 1))
C02 = inner(Dx(u[1], 1), Dx(v[0], 0))
C03 = inner(Dx(u[1], 1), Dx(v[1], 1))
C = [C00, C01, C02, C03]

# Compute BlockMatrix
M = BlockMatrix(A + C)

# Compute right hand side
# Use sympy to compute a rhs
x, y = symbols("x,y")
fe = (sin(x**2)*(1-y)*(1+y), y**2*(cos(x-y**2)))
fj = Array(V, buffer=fe)

b = inner(v, fj)

# Compute solution
u_hat = M.solve(b)
u_sol = u_hat.backward()

[mikaem]

Hi Thilo

Thanks a lot for reporting. This is a bug. The form represents a scalar and it should not go into this if test. I know this used to work, so this only goes to show that we need better tests! I'll look into it right away.

[mikaem]

If you want to fix this locally, then the following added just above mentioned line should do:

if test.expr_rank() == 0:
    recursive = False

[sebglane]

Hi Mikael,

Thilo and I realized that you added the lines above in commit 955a326ba71e65d75968cbfda2ad0706a0cca8c1. The code which is shown above is working properly now. Thank you for fixing this.

However, we think that there is deeper a problem somewhere. If we modify the boundary conditions in our example from homogeneous Dirichlet bcs to inhomogeneous ones, an error occurs if the linear is solved. We hope that example below show the problem. We are not sure if we should close this issue and open a new one because the problem does not seem to be related to the div-div-term. The problem also occurs if only the grad-grad is considered.

Best wishes, Sebastian and Thilo

from sympy import symbols, cos, sin
from mpi4py import MPI
from shenfun import inner, grad, div, TestFunction, TrialFunction, Array,\
    Basis, TensorProductSpace, VectorTensorProductSpace, BlockMatrix
from mpi4py_fft.pencil import Subcomm

comm = MPI.COMM_WORLD

family = 'legendre'

# Size of discretization
n = 40
N = [n, n+1]

# Bases and spaces
BX_hom = Basis(N[0], family=family, bc=(0, 0))
BX_inhom = Basis(N[0], family=family, bc=(0, 0.1))
BY_hom = Basis(N[1], family=family, bc=(0, 0))

# Function spaces
subcomms = Subcomm(MPI.COMM_WORLD, [0, 1])
T_hom = TensorProductSpace(subcomms, (BX_hom, BY_hom))
T_inhom = TensorProductSpace(subcomms, (BX_inhom, BY_hom))

V_hom = VectorTensorProductSpace([T_hom, T_hom])
V_inhom = VectorTensorProductSpace([T_inhom, T_inhom])
V_mixed = VectorTensorProductSpace([T_hom, T_inhom])

# Choose the BCs you want
# working
V = V_hom

# not working
# V = V_inhom
# V = V_mixed

# Functions
u = TrialFunction(V)
v = TestFunction(V)

# Compute matrices
# Assemble standard term
A = inner(grad(u), grad(v))

# Assemble (div(u), div(v))-term
# non-working
B = inner(div(u), div(v)) # raises an error

# Compute BlockMatrix
M = BlockMatrix(A + B)

# Compute right hand side
# Use sympy to compute a rhs
x, y = symbols("x,y")
fe = (sin(x**2)*(1-y)*(1+y), y**2*(cos(x-y**2)))
fj = Array(V, buffer=fe)

b = inner(v, fj)

# Compute solution
u_hat = M.solve(b) # This line fails because inconsistent shapes
u_sol = u_hat.backward()

[mikaem]

Hi

Non-homogeneous bcs are tricky. Please have a close look at, e.g., the driven cavity demo. There are matrices in A and B that should go on the right hand side. Please also look at the Rayleigh Bénard demo. Best of luck😀 Mikael