PDE in divergent form

[erkara]

Hello. I am trying to solve the following PDEs representing volume and force balance in a porous medium.

$$ \nabla \cdot \left( \frac{\partial \mathbf{u}}{\partial t} - \mathbf{K}\nabla P \right) = f(\mathbf{x},t) $$

$$ \nabla \cdot \left( \mathbf{\sigma} - P\mathbf{I} \right) = 0 $$

where $\mathbf{u}$ is the vector valued displacement and and $P$ is the pressure while $\mathbf{K}$ is permeability tensor and $\mathbf{\sigma}= \mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) + \lambda (\nabla \cdot \mathbf{u}) \mathbf{I}$ is the solid stress term. It is given that $u=0$ and $P=400$ on the boundary $\partial\Omega$. Using scalar test functions $\phi$ and $\psi$ vanishing on the boundary, we obtain the following weak forms:

$$F_u:=- \int_{\Omega} \Big[ \nabla \phi \cdot \left( \frac{\partial \mathbf{u}}{\partial t} +\mathbf{K} \nabla P \right) + \phi f(\mathbf{x},t) \Big]d\Omega = 0$$

$$F_p:= \int_{\Omega} \Big[-\nabla \psi : \mathbf{\sigma} + (\nabla \cdot \psi) P \Big] d\Omega = 0$$

I am trying to set up this problem on Firedrake with simplest case possible with $K = \mathbf{I}$ and $f=1$. Assuming the weak form is correct and I am not missing anything, I set up the problem as follows. I am new to Firedrake and Fenics, my understanding is that solution vector $\mathbf{u}$ and the scalar test function $\phi$ should use the same \textit{VectorFunctionSpace} object but this naturally gives me "Shapes do not match:" error in the variational form from the line I define the variational form $F_u$a as well as $F_p$. Should I somehow pull these test functions from something like Q = FunctionSpace(mesh,"CG",degree=1)? Here is my code. I highly appreciate any kind of insight:

from firedrake import *


def main():
    # some simulations parameters and maybe more
    T = 1
    dt = T/60

    # Create mesh and define mixed function space
    nx, ny = 64, 64
    mesh = RectangleMesh(nx, ny, 5, 5, quadrilateral=True)
    V = VectorFunctionSpace(mesh,"CG",degree=2)
    Q = FunctionSpace(mesh,"CG",degree=1)
    Z = MixedFunctionSpace([V, Q])  # Mixed function space

    v,q = TestFunctions(Z)

    #current and previous solutions for u,p and c
    U = Function(Z)
    U_n = Function(Z)

    #init set-up
    t = Constant(0.0)
    x = SpatialCoordinate(mesh)

    #boundary conditions
    bc_u = DirichletBC(Z.sub(0), Constant((0, 0)), "on_boundary")
    bc_p = DirichletBC(Z.sub(1), Constant(400), "on_boundary")


    bcs = [bc_u, bc_p]

    # apply initial conditions
    u_ic = Function(Z.sub(0)).interpolate( Constant((0, 0))  )
    p_ic = Function(Z.sub(1)).interpolate(Constant(50))
    U_n.sub(0).assign(u_ic)
    U_n.sub(1).assign(p_ic)

    #Split mixed functions
    u, p = split(U)
    u_n, p_n = split(U_n)

    # Volume balance--> incompetible function space error
    f = Constant(1.)
    K = Constant(((1, 0), (0, 1)))
    u_diff = (u - u_n) / dt

    F_u = -(inner(u_diff, v) + inner(grad(v), K * grad(p)) + dot(v, f))*dx

    # displacement--> incompetible function space error
    E = Constant(5000)
    nu = Constant(0.4)
    lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
    mu = E / (2 * (1 + nu))
    F_p = (inner(grad(q), mu * (grad(u) + grad(u).T)) + div(q) * lambda_ * div(u) - div(q) * p) * dx

    #nonlinear solve stuff etc

if __name__ == "__main__":
    main()

[wence-]

To make the shapes match up, it looks like you have a Petrov-Galerkin, rather than Galerkin method. So you want to find a vector-valued $u$ and scalar-valued $P$. But your test function $\phi$ must be scalar-valued, while $\psi$ must be vector-valued. Firedrake doesn't have great support for Dirichlet conditions in Petrov-Galerkin methods, so I am not sure this would work, but something that will at least allow you define the forms would be:

# Note: I have no idea what the right sets of spaces are
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)

Z = V*Q
z = Function(Z)

u, p = split(z)

# This is probably the wrong pair, but the shapes are right
W = Q*V
phi, psi = TestFunctions(W)

I = Identity(2)
sigma = mu*(grad(u) + grad(u).T) + lambda_ * div(u) * I
F = (
    inner(grad(phi), u_diff) 
    + phi * f
    - inner(grad(psi), sigma)
    + div(psi)*p
)*dx

bcs = [DirichletBC(Z.sub(0), 0, "on_boundary"), DirichletBC(Z.sub(1), 400, "on_boundary")]

solve(F == 0, z, bcs=bcs)

I think the bcs won't work properly since Firedrake won't correctly identify the columns in the operator that match up with the boundary rows. See #197

[dham]

Your first issue is that you are assuming K *grad(p) does a matrix-vector multiplication, which it doesn't. The line should instead be:

i, j = indices(2)
F_u = -(inner(u_diff, v) + grad(v)[i] * K[i,j] * grad(p)[j] + dot(v, f))*dx

This uses Einstein summation to represent the two contractions in your maths.

The issue with F_p is a problem in your maths which is being continued through in the code. If q/$\psi$ is a scalar test function then grad(q) is a vector. However, u is a vector so grad(u) + grad(u).T is a dimension 2 tensor and the inner product you have written is not well-defined. This issue is already in the maths so you need to look back to your source and work out what is going on.

Incidentally, another issue in your maths is that v/$\phi$ is not a scalar test function but is vector-valued (because u is a vector and hence the corresponding test space is vector valued). This one isn't causing you the same issue in the code because you do actually treat v as a vector despite what is written above.

[erkara]

Thanks for your reply @dham. You are right that $\psi$ must be a vector test function. I guess I could do a similar double contraction there. However, I still cannot compile F_u as you suggested as I get "Component and shape length don't match"

As pointed out by @wence-, what is the right functions spaces to use here? I guess this is more of a FEM question than Firedrake.

[dham]

I think the error is still in your weak form. You need to carefully think through all the shapes there so that your integrands are all scalar. The Firedrake errors you are seeing are just the result of that.