Why do the shapes below not match?

Question

Hello. I am trying to implement a variational problem involving pressure and one-dimensional displacement. I get the following error from the code:

Shapes do not match: <Grad id=139643460028464> and <Indexed id=139643460375568>.

The error occurs at the lines where the variational problem is defined ("F = p_v_1dx + ...").

Here is the relevant code:

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# Parameters
c = 1.25
l = 5.0 
P_surf = 0.0
sigma_0 = 3.0 
gamma = 0.65 
P0 = sigma_0*gamma 
Kv = 10 
cm = 5
# Derivative of pressure at the bottom surface.
der_p_bottom = 0.0
# Derivative of displacement at the top surface.
der_w_top = sigma_0/Kv

# The initial time.
t = 0.0
# The final time
T = 2.0
# The number of time steps to take.
num_steps = 40
# The step size
dt = (T-t)/num_steps


# Create mesh and define function space.
mesh = IntervalMesh(100,0,l)
P1 = FiniteElement('P', interval, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define boundary conditions for pressure.
tol = 1E-14

def upper_boundary(x,on_boundary):
    return on_boundary and abs(x[0])<tol

# An expression for the derivative of P at the bottom surface.
der_p_D = Expression('der_p_bottom*x[0]/l',degree = 1, \
    der_p_bottom=der_p_bottom, l=l)

bc_p = DirichletBC(V.sub(0),Constant(P_surf),upper_boundary)

# Define boundary conditions for displacements.
def bottom_boundary(x,on_boundary):
    return on_boundary and abs(x[0]-l)<tol

der_w_D = Constant(der_w_top)

bc_w = DirichletBC(V.sub(1),Constant(0.0),bottom_boundary)

# Aggregate the Dirichlet boundary conditions in a list.
bc = [bc_p, bc_w]

# Define the test function space.
v_1, v_2 = TestFunctions(V)

# Define functions for solutions at previous and current time steps.
# Previous time step.
# Define the initial conditions.
init_cond = Expression(('P0','sigma_0/Kv*(l-x[0])'), degree = 1,
                       sigma_0=sigma_0, Kv=Kv, l=l,P0=P0)
u_n = interpolate(init_cond,V)
p_n,w_n = u_n.split()            
# Current time step.
u = Function(V)
p_, w_ = split(u)

# Define the variational problem.
F = p_*v_1*dx + c*dt*inner(grad(p_),grad(v_1))*dx - c*dt*der_p_D*v_1*ds \
    - p_n*v_1*dx \
    + inner(grad(w_),grad(v_2))*dx + cm*inner(grad(p_),v_2)*dx-der_w_D*v_2*ds

# Change the variational problem into a bilinear form.
a,L = lhs(F),rhs(F)

# Now, declare the current time step's variable as a function.
p_ = Function(V)

# Create the progress bar.
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Open the file to store results in VTK format.
vtkfile_w = File('Terzaghi_results_FI/pressure.pvd')

# Form and solve the linear system.
for n in range(num_steps):

    # Update current time.
    t += dt

    # Solve the variational problem.
    solve(a == L, u, bc)

    # Extract the solution.
    p_,w_ = u.split()
    p_comp = p_.compute_vertex_values(mesh)
    w_comp = w_.compute_vertex_values(mesh)

    # Write the results to a VTK file for visualization in Paraview.
    vtkfile_w << (p_,t)
    vtkfile_w << (w_,t)

    # Update the previous solution.
    u_n.assign(u)

    # Update the progress bar.
    progress.update(t/T)

Any thoughts on why this may be occurring?

Answer 1

It is easy test out the individual terms: It is the second last term, which is causing the issue.

I suggest replacing

 cm*inner(grad(p_),v_2)*dx

with

 cm*p_.dx(0)*v_2*dx

Comments

KristianE, thank you for your response. I appreciate it. The error no longer occurs at the above line. Two follow up questions:

1. How can one check the shapes of the individual terms? I am using
   Spyder as my editor, and the variable window does not show functions
   such as p_ and w_.
2. Regarding your expression above, is p_.dx(0) equivalent to grad(p_)?
   What is the .dx property?

---

.dx(0) is the first (x) spatial derivative. An alternative solution, is to use a VectorFiniteElement for the velocity. This is probably the better solution, since you would need to that anyway for dimensions higher than 1.

---

Thank you for your answer. The code now works.