I am trying to create a fully-implicit solver for a system of two PDEs. I am using a mixed element consisting of a vector element for displacement and a scalar-valued element for pressure. I have been looking at the examples in the FEniCS tutorial (the Navier-Stokes example and the advection-reaction diffusion equations), but when I run my code, I get the following error:
Can only integrate scalar expressions. The integrand is a tensor
expression with value shape (2,) and free indices with labels ().
The relevant code:
from fenics import *
import numpy as np
# Parameters
c = 1.25
l = 5.0
P_surf = 0.0
sigma_0 = 3.0
gamma = 0.65
P0 = sigma_0*gamma
Kv = 10.0
cm = 5.0
# Derivative of pressure at the bottom surface.
der_p_bottom = 0.0
# Derivative of displacement at the top surface.
der_w_top = cm*gamma*sigma_0
# 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
# Width
width = 3.0
# Thickness
thickness = 2.0
# Create mesh and define function space.
# Mesh for the 2-D case.
mesh = RectangleMesh(Point(0.0,0.0),Point(l,width),100,10)
# Use a vector-valued element for displacement.
P1 = VectorElement('P', triangle, 1)
# Use a scalar-valued element for pressure.
P2 = FiniteElement('P', triangle, 1)
# Create a mixed element to solve the system using the fully-implicit method.
element = MixedElement([P1, P2])
# Create the function space corresponding to this mixed element.
V = FunctionSpace(mesh,element)
# Create the test functions.
v_1, v_2 = TestFunctions(V)
# Create the trial functions.
ans = Function(V)
w_, p_ = split(ans)
# Define boundary conditions for pressure.
upper_boundary = 'on_boundary && near(x[0],0)'
# An expression for the derivative of P at the bottom surface.
der_p_D = Expression(('near(x[0],l) ? der_p_bottom : 0.0','0.0'),degree = 1, \
der_p_bottom=der_p_bottom, l=l)
# Set the pressure at the porous medium's top to 0 (second subsystem of V).
bc_p = DirichletBC(V.sub(1),Constant(P_surf),upper_boundary)
# Define boundary conditions for displacements.
bottom_boundary = 'on_boundary && near(x[0],%s)'%(l)
der_w_D = Expression(('near(x[0],0.0) ? der_w_top : 0.0','0.0'),degree=0, \
der_w_top=der_w_top)
# Set the displacement at the porous medium's bottom to 0 (first dimension of
# the first subsystem of the FunctionSpace V).
bc_w = DirichletBC(V.sub(0).sub(0),Constant(0.0),bottom_boundary)
# Define functions for solutions at previous and current time steps.
# Previous time step.
# Define the initial conditions.
init_cond = Constant(P0)
p_n = interpolate(init_cond,V.sub(1).collapse())
# Define the variational problem.
# Equation 1
Eq1 = p_*v_1*dx - p_n*v_1*dx \
+c*dt*inner(grad(p_),grad(v_1))*dx \
-c*dt*der_p_D*v_1*ds
# Equation 2
Eq2 = inner(grad(w_),grad(v_2))*dx \
-der_w_D*v_2*ds \
+cm*p_.dx*v_2*dx
# Add the equations to make a monolithic scheme.
F = Eq1 + Eq2
# Create the progress bar.
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
# Form and solve the linear system.
for n in range(num_steps):
# Update current time.
t += dt
# Solve the variational problem.
solve(F == 0, ans, bc)
# Check on the solution.
w_,p_ = ans.split()
p_comp = p_.compute_vertex_values(mesh)
w_comp = w_.compute_vertex_values(mesh)
# Update the previous solution.
p_n.assign(p_)
# Update the progress bar.
progress.update(t/T)
The error occurs at the line "-cdtder_p_Dv_1ds" in defining Eq1. Have the test functions been made properly? If yes, what could be the problem in this case?
