I have solved the Helmholtz equation with complex “wave number” in dolfin without problems by splitting it in its real and imaginary parts and using mixed function spaces.
The problem along with the code is already shown here, but I copy here the code, just in case:
#imports
from dolfin import RectangleMesh, MeshFunction, SubDomain, Measure, FunctionSpace, MixedFunctionSpace, Constant
from dolfin import TrialFunctions, TestFunctions, inner, grad, dx, Function, solve, plot, Point, near, Expression
from dolfin import assemble
#geometrical constants and mesh generation
Lx = 0.01
aspect_ratio = 100
Ly = aspect_ratio*Lx
Nx = 5
Ny = aspect_ratio*Nx
mesh = RectangleMesh( Point(0,0), Point(Lx,Ly), Nx, Ny )
dim_edge = 1
def_value = 0
boundary_markers = MeshFunction("size_t", mesh, dim_edge, def_value)
tol = 1.0E-14
front = 1
back = 2
class Gamma_front(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)
GF = Gamma_front()
GF.mark(boundary_markers, front)
class Gamma_back(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], Lx, tol)
GB = Gamma_back()
GB.mark(boundary_markers, back)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
#Physical constants
kappa = 200
h_conv = 15
epsilon = 0.8
sigma = 5.67E-8
T_air = 298
rho = 2700
c= 900
omega = 0.8
alpha = 0.4
I = 600*(1 + 0.0j)
sx = -15*Lx
sy = 0.5*Ly
Q_R = Expression("alpha*I_R/(2*pi)*(x[0]-sx)/((x[0]-sx)*(x[0]-sx) + (x[1]-sy)*(x[1]-sy))", degree=2, alpha = alpha, I_R = I.real, sx=sx, sy=sy)
Q_I = Expression("alpha*I_I/(2*pi)*(x[0]-sx)/((x[0]-sx)*(x[0]-sx) + (x[1]-sy)*(x[1]-sy))", degree=2, alpha = alpha, I_I = I.imag, sx=sx, sy=sy)
# Funtion spaces
P2 = FunctionSpace(mesh, "Lagrange", 1 )
H = MixedFunctionSpace( P2, P2 )
Phi_R, Phi_I = TrialFunctions(H)
Psi_R, Psi_I = TestFunctions(H)
# variational forms
a = kappa*( inner( grad(Phi_R), grad(Psi_R)) + inner( grad(Phi_I), grad(Psi_I) ))*dx \
- rho*omega*c*( Phi_R*Psi_I - Phi_I*Psi_R)*dx \
+ h_eff*( Phi_R*Psi_R + Phi_I*Psi_I)*ds(front) \
+ h_eff*( Phi_R*Psi_R + Phi_I*Psi_I)*ds(back)
L = ( Q_R*Psi_R + Q_I*Psi_I )*ds(front)
#problem solution
T = Function(H)
solve(a==L, T)
T_R, T_I = T.split()
However, now I was interested in solving similar problems with different right hand sides, so I tried to compute the assembled matrix of the bilinear form using:
from dolfin import assemble
A = assemble(a)
b = assemble(L)
and I get the message error: ArityMismatch: Adding expressions with non-matching form arguments ('v_0^0', 'v_1^1') vs ('v_0^1', 'v_1^0').
After performing a look on the internet I saw that this can be caused by the forms not having all the trial and test functions of the mixed space, but I don’t understand how to fix it.
Thank you very much