Hello, I’m trying to solve the Helmholtz equation but with a complex coefficient (and hence solving it on complex functions). I obtained the variational formulation in terms of u_r and u_i, the real and imaginary parts of my solution, but when implementing it on FEniCS I get an exception with error message:
UFLException: Invalid ranks 1 and 1 in product.
Here is what I think is the MWE:
# parameters
lx = 0.01
ly = 1
k = 200
h = 15
alpha = 0.5
rho = 2700
c = 890
Nx = 5
Ny = 100
from dolfin import *
# Create mesh and define function space
mesh = RectangleMesh( Point(0, 0), Point(lx, ly), Nx, Ny, diagonal= 'right' )
P2 = VectorElement('P', triangle, 1)
Element = P2*P2
W = FunctionSpace(mesh, Element)
# Define boundary conditions
q_r = Expression('1000*(0.2)/( (0.2)*(0.2) + (x[1]-0.5)*(x[1]-0.5))', degree=2)
q_i = Expression('1000*(0.2)/( (0.2)*(0.2) + (x[1]-0.5)*(x[1]-0.5))', degree=2)
tol = 1E-14
class BoundarySides(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ( near(x[1], 0, tol) or near(x[1], ly, tol) )
class BoundaryFront(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)
class BoundaryBack(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], lx, tol)
sides = BoundarySides()
front = BoundaryFront()
back = BoundaryBack()
boundary_markers = MeshFunction("size_t", mesh, 1, 0)
#labeling the markers with different numbers
sides.mark(boundary_markers, 0)
front.mark(boundary_markers, 1)
back.mark(boundary_markers, 2)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Define variational problem
#U = Function(W)
#u_r, u_i = U.split()
(u_r, u_i) = TrialFunctions(W)
(v_r, v_i) = TestFunctions(W)
a = k*inner(nabla_grad(u_r), nabla_grad(v_r))*dx + k*inner(nabla_grad(u_i), nabla_grad(v_i))*dx \
+ rho*c*u_i*v_r*dx - rho*c*u_r*v_i*dx \
+ h*u_r*v_r*ds(1) + h*u_r*v_r*ds(2) + h*u_i*v_i*ds(1) + h*u_i*v_i*ds(2)
L = alpha*q_r*v_r*ds(1) + alpha*q_i*v_i*ds(1)
#F = a - L
# Compute solution
T_num = Function(W)
T_r, T_i = T_num.split()
solve(a == L, T_num)
#solve(F==0,T_num)
As you can see I have tested also solving it as a nonlinear problem but I got the same error.
The error I get points to the second line of the definition of a
that is, when the scalar multiplication of u_r
and v_r
takes place.
Thank you very much for you help