Problem with bc.apply()

1

Hi all,
I have tried to follow the Naiver-stoker example in the tutorial and solve two coupled heat transfer and electrostatic PDEs separately. However, the following code gives the following error:

File “thermoelectric.py”, line 98, in
bc_potential.apply(A2)
TypeError: apply(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericVector) -> None
2. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericMatrix) -> None
3. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericMatrix, arg1: dolfin::GenericVector) -> None
4. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericVector, arg1: dolfin::GenericVector) -> None
5. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericMatrix, arg1: dolfin::GenericVector, arg2: dolfin::GenericVector) -> None

The problem is the line where I apply the bc, bc_potential.apply(A2) and also the other one, give anyone please explain why the code is complaining about this line? The code is attached below, many thanks for your help.

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


box = Box(Point(-0.3,-0.15,-0.014),Point(0.3,0.15,0.014)) ## sample dimensions in cm
mesh = generate_mesh(box,100)
P1 = FiniteElement('P',tetrahedron,1)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh, element)

T = 0.001
num_steps = 10
dt = T/num_steps
a = 0.6 ## sample length in cm
b = 0.3 ## sample width in cm
c = 0.028 ## sample thickness in cm
L = 100.0e-4 ## wire lenght in cm
r = 5.0e-4 ## wire radius in cm
R = (L/(r*r*np.pi))*(2.65e-4) ## wire resistance in ohm
I = 200e-3 ## current in A
contact = 4e3 ## contact resistance in ohm
W = I*I*(R+contact) ## heat injected via the wire
T_ext = 30.0 ##environment temperature in K
J_ext = I/(r*r*np.pi) ## current density in wire in A/cm2
tol = 0.05

# Define boundary condition
u_D_heat = Expression('T_ext', degree=2,T_ext=T_ext)
u_D_potential = Expression('near(x[0],-0.25,tol) ? 100 : near(x[0],0.25,tol) ? -100 : 0',degree=2,tol=tol)
def boundary_D_heat(x, on_boundary):
    return False
def boundary_D_potential(x, on_boundary):
    return False

bc_heat = DirichletBC(V.sub(0), u_D_heat, boundary_D_heat)
bc_potential = DirichletBC(V.sub(1), u_D_potential, boundary_D_potential)
#bc = [bc_heat, bc_potential]

# Define initial condition
u_0 = Expression(('30.0','0.0'),degree=2)
u_n = interpolate(u_0,V)
u_n1, u_n2 = u_n.split()

# Define parameter functions
K = Expression('temp < 20 ? 0.0011*pow(temp,4)-0.0567*pow(temp,3)+0.978*pow(temp,2)-3.0736*temp+3.1126 : 0.00000002*pow(temp,4)-0.00002*pow(temp,3)+0.0064*pow(temp,2)-1.0824*temp+69.056',degree=4,temp=u_n1)
#K = Constant(0.13) # w/cm/K
cp = Constant(0.2) #J/g/K
rou = Constant(2.329) #g/cm3
sigma = Expression('temp < 37.0 ? exp(0.4258*temp-23.313):exp(0.0105*temp-7.6668)',degree=2,temp=u_n1) #conductivity
h = 0.002 # w/cm2/K
k = Constant(dt)
# Define variational problem
u = Function(V)
u_1, u_2 = split(u)
v_1, v_2 = TestFunctions(V)

g_heat = Expression('(near(x[0],-0.25*a,tol) and near(x[1],0.0,tol) and near(x[2],c/2.0,tol/5.0)) or (near(x[0],0.25*a,tol) and near(x[1],0.0,tol) and near(x[2],c/2.0,tol/5.0)) ? W : h*(T_ext-temp)',degree=2, tol=tol,a=a,c=c,W=W, h=h,T_ext=T_ext,temp=u_n1)
g_potential = Expression('near(x[0],-0.25*a,tol) and near(x[1],0.0,tol) and near(x[2],c/2.0,tol/5.0) ? -1.0*J_ext : near(x[0],0.25*a,tol) and near(x[1],0.0,tol) and near(x[2],c/2.0,tol/5.0) ? 1.0*J_ext : 0',degree=2,tol=tol,a=a,c=c,J_ext=J_ext)
#g_potential = Constant(0.0)
f_heat = sigma*dot(grad(u_2),grad(u_2))

#electrostatic problem
F2 = sigma*dot(grad(u_2),grad(v_2))*dx - g_potential*v_2*ds
a2 = lhs(F2)
L2 = rhs(F2)

#heat transfer problem
F1 = k*K*dot(grad(u_1), grad(v_1))*dx + rou*cp*u_1*v_1*dx - k*f_heat*v_1*dx - rou*cp*u_n1*v_1*dx - k*g_heat*v_1*ds
a1 = lhs(F1)
L1 = rhs(F1)


A1 = assemble(a1)
A2 = assemble(a2)

bc_potential.apply(A2)
bc_heat.apply(A1)


# Time-stepping
t = 0.0
vtkfile_u_1 = File('heat_test/solution_u_1.pvd')
vtkfile_u_2 = File('heat_test/solution_u_2.pvd')

for n in range(num_steps):
    t = t + dt

    b2 = assemble(L2)
    bc_potential.apply(b2)
    solve(A2, u_2.vector(), b2)

    b1 = assemble(L1)
    bc_heat.apply(b1)
    solve(A1, u_1.vector(), b1)

    u_n.assign(u)
    _u_1, _u_2 = u.split()
    vtkfile_u_1 << (_u_1,t)
    vtkfile_u_2 << (_u_2,t)
#c = plot(u, mode='color')
#plt.colorbar(c)
#plt.show()
2

You should not use mixed functionspaces when trying to solve problems separately. See the navier-stokes demo which uses a splitting scheme.
Here there is usage of TrialFunction and TestFunction for each PDE and the solution for each PDE is given a specific variable (u0, u1, p1).

The specific error you got was because you tried to use assemble a variational form which did not have a trial function in it, yielding lhs(F2)=0.

3

Thanks for your reply, dokken. I did not use the TrialFunction because I learned that for nonlinear PDE, one should not use TrialFunction instead, Function should be used, maybe I did not really understand the difference between TrialFunction and Function. So following your suggestion, I use Trial function for the definition of the variational forms and then define the Function, the problem is solved. Thanks.