Shape mismatch error

Question

Hi,
I'm trying to solve the following system of equations on a 1D interval mesh:

$$\frac{\partial A}{\partial t} + \nabla \left( A u \right) = 0$$
$$\frac{\partial u}{\partial t} + \frac{\beta}{2 \rho \sqrt{A}} \nabla A + u \nabla u = - \frac{8 \pi \mu}{\rho} \frac{u}{A}$$

The equations describe the time evolution of the cross-sectional area of an artery and blood velocity u within that artery. I've gathered from the tutorial examples that a useful approach for solving a highly coupled system like this is the use of two test functions v and q. I use v as the test function for the first equation and q as the test function for the second equation and then add both equations so that I am only left with a single equation (see code below). As I implement a(u, A, v, q) I run into a shape mismatch error, which occurs on the term dt*beta/(2*rho)*inner(a**-0.5*grad(a), q)*dx. How do I resolve this error?

# Set up simulation parameters
l = 20.0
a0 = 1.0
beta = 22967.4
c0 = 357.14
mu = 3.5e-3
rho = 1060.0

nx = 100
dx = l/nx
dt = 0.5*dx/c0
T = 0.06
nt = int(T/dt)
nu = 3.5e-3


# Create 1d line mesh
mesh = IntervalMesh(nx, 0.0, l)

# Define function spaces (P2-P1)
V = FunctionSpace(mesh, "DG", 1)
Q = FunctionSpace(mesh, "DG", 1)
W = MixedFunctionSpace([V, Q])

# Define trial and test functions for velocity and area
(u, a) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# Define boundary conditions
p_in = 1.0
p_out = 0.0
a_in = A0 / (1 - p_in/(2*rho*c0**2))**2
a_out = A0 / (1 - p_out/(2*rho*c0**2))**2
L = Constant(l)
inflow  = DirichletBC(W.sub(1), Constant(a_in), 
                      "on_boundary && near(x[0], 0.0)")
outflow = DirichletBC(W.sub(1), Constant(a_out), 
                      "on_boundary && near(x[0], 1.0)")
bcs = [inflow, outflow]

u_0 = Function(W)
a_0 = Function(W)

a = inner(a, v)*dx + dt*inner(u*grad(a), v)*dx + dt*inner(a*grad(u), v)*dx +\
    inner(u, q)*dx + dt*beta/(2*rho)*inner(a**-0.5*grad(a), q)*dx +\
    dt*inner(u*grad(u), q)*dx
L = inner(a_0, v)*dx + inner(u_0, q)*dx -\
    8*np.pi*mu*dt/rho*inner(u/a, q)*dx
A = assemble(a)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True

t = dt
while t < T + DOLFIN_EPS:
    b = assemble(L)
    [bc.apply(A, b) for bc in bcs]
    solve(A, (u_0.vector(), a_0.vector()), b, "bicgstab", "default")

    plot(u_0, title="velocity", rescale=True)
    plot(a_0, title="area", rescale=True)    

    t += dt
    print "t=", t

# Hold plot
interactive()

Comments

Can you explain the term: a**-0.5*grad(a)?

---

Hi,
sorry, I thought it was obvious from the equations and the code. It's the variational form of the term $$\frac{\beta}{2 \rho} \frac{\nabla A}{\sqrt{A}}$$ in the second equation, where the TrialFunction a in the code is the variable A in the equation.

Answer 1

Since u and q are vectors, you need a vector functionspace and I think in this problem you don't need to use discontinuous elements so: use Q = VectorFunctionSpace(mesh, "CG", 1).

The weak formulation also has some issues. The weak form of the first equation is:

 inner(a, v)*dx + dt*inner(grad(a*u),v)*dx 

or if you want to integrate by parts then:

 inner(a, v)*dx + dt*inner(a*u,grad(v))*dx - (boundary terms)

Comments

Thank you for your suggestions! I have made the changes you suggested, ie have corrected the weak form of the first equation to

inner(a, v)*dx + dt*inner(grad(a*u), v)*dx

and have changed my function spaces

V = VectorFunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, Q])

But I still get the same shape mismatch error on the same line.

---

Are you sure, that the problem is in same line? If I run your code I get an error related to the boundary conditions.

---

I can reproduce the boundary condition error if I define both V and Q as a VectorFunctionSpace, but thinking about it, shouldn't they both be FunctionSpace as I'm working in 1D.

---

In my understanding velocity is a vector even in 1D (it has magnitude and direction).

Answer 2

You are redefining a from being a trial function to a trial test function product. Furthermore, last I checked, one should not use trial functions for solving non-linear problems.