UFL: unsupported operand type(s) for +: 'Sum' and 'Form'

Question

I'm trying to write a solver to recreate the results of the following 1D polymer curing problem:
https://www.comsol.com/blogs/modeling-the-thermal-curing-process/

But I get the error:
TypeError: unsupported operand type(s) for +: 'Sum' and 'Form'
on the variational statement line

What does the error mean and how could I fix it?

The full code is shown below:

from dolfin import *

# Constants (all in SI units)
sf = 1.0e0 #length scale factor (m to m)

rho = 1200.0/sf/sf/sf #density
C_p = 1000.0*sf*sf    #specific heat
k   = 0.2/sf          #thermal conductivity
H_r = 5.0e5*sf*sf     #total heat of reaction
k_A = 20.0e3          #frequency factor
E_a = 50.0e3*sf*sf    #activation energy
R   = 8.314*sf*sf     #universal gas constant
n   = 1.4             #order of reaction
l   = 5.0e-3*sf       #length of 1D interval
numElems = 100        #number of finite elements
Q_heat   = 10.0e3     #10 kW/m^2 for one m^2 which is 10 kW of heating power
dt = 1.0     


# Define function space for system of concentrations
mesh = IntervalMesh(numElems, 0.0,l)
P1 = FiniteElement('Lagrange', interval, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)


# Define test functions
v, v_2 = TestFunctions(V)

# Define functions for velocity and concentrations
u   = Function(V)
u_n = Function(V)

# Split system functions to access components
u, alpha     = split(u)
u_1, alpha_1 = split(u_n)


# Initial conditions
u_0 = Expression(('0.0', '20.0+273.15'), degree=1)
u_n = project(u_0, V)


# Boundary marker
class HeatFlux(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[0] > (l-DOLFIN_EPS))

# Initialize sub-domain instances
right = HeatFlux()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
right.mark(boundaries, 1)
dss = ds[boundaries]


# Define source term
f = Expression('0.0',degree=1)

# Define expressions used in variational forms
Dt  = Constant(dt)
k_A = Constant(k_A)
k2  = Constant(E_a/R)
n   = Constant(n)
rho = Constant(rho)
C_p = Constant(C_p)
k   = Constant(k)
H_r = Constant(H_r)


F = alpha - alpha_1 + Dt*k_A*exp(-k2/u)*(1-alpha)**n + \
    alpha - alpha_1 - f*Dt + \
    rho*C_p*u*v*dx + Dt*k*inner(nabla_grad(u), nabla_grad(v))*dx - (rho*C_p*u_1 - Dt*rho*H_r*f)*v*dx + Dt*g*v*dss(1)


T = 2.0 # Total simulation time: 8 minutes
t = 0.0
while t <= T:
    print 'time = %f' % t

    t += dt

    # Solve variational problem for time step
    solve(F == 0, u)

    # Update previous solution
    u_n.assign(u)


plot(u, title="Temperature")
interactive()

Answer 1

There are components of your residual which are not correctly formulated:

alpha - alpha_1 + Dt*k_A*exp(-k2/u)*(1-alpha)**n + \
alpha - alpha_1 - f*Dt + \

These need to be tested against the function v and integrated over the volume dx or boundary ds correctly.

Comments

Thank you for your answer Nate. I updated the program but I still get the error:

Error: Unable to solve variational problem.
*** Reason: Expecting second argument to be a Function.

This is in the solver. I updated the program and added another trial function based on alpha. The code is shown below. I was able to get this working by solving two different systems. I think this method will be faster but I would like to get this mixed formulation method working to see it.

The method using two different solves is:

while t <= T:
    print 'time = %f' % t
    b = assemble(L, tensor=b)
    solve(A, u.vector(), b)

    t += dt

    u_1.assign(u)

    G = (alpha - alpha_1 - dt*k_A*exp(-k2/u)*(1-alpha)**n)*v_2*dx
    solve(G == 0, alpha)

    f = (alpha-alpha_1)/dt 
    L = (rho*C_p*u_1 - dt*rho*H_r*f)*v*dx + dt*g*v*ds(1)

    alpha_1.assign(alpha) 

This method should be faster because: the heat equation is linear in space (within the time loop) and can be solved with a linear solver and the matrix A for the heat equation is constant so it needs to be assembled once, then the nonlinear equation with alpha can be solved on its own as a nonlinear problem.

The issue here is trying to implement the statment f = (alpha-alpha_1)/dt in the mixed formulation. It has terms from the original heat equation and the alpha equation but no variables of its own. Do I need a third trial function for this equation?

The updated mixed program (with the trial function for alpha) is shown below:

from dolfin import *

# Constants (all in SI units)
sf = 1.0e0 #length scale factor (m to m)

rho = 1200.0/sf/sf/sf #density
C_p = 1000.0*sf*sf    #specific heat
k   = 0.2/sf          #thermal conductivity
H_r = 5.0e5*sf*sf     #total heat of reaction
k_A = 20.0e3          #frequency factor
E_a = 50.0e3*sf*sf    #activation energy
R   = 8.314*sf*sf     #universal gas constant
n   = 1.4             #order of reaction
l   = 5.0e-3*sf       #length of 1D interval
numElems = 100        #number of finite elements
Q_heat   = 10.0e3     #10 kW/m^2 for one m^2 which is 10 kW of heating power
dt = 1.0     


# Define function space for system of concentrations
mesh = IntervalMesh(numElems, 0.0,l)
P1 = FiniteElement('Lagrange', interval, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)


# Define test functions
v, v_2 = TestFunctions(V)

# Define functions for velocity and concentrations
u   = Function(V)
u_n = Function(V)

# Split system functions to access components
u, alpha     = split(u)
u_1, alpha_1 = split(u_n)


# Initial conditions
u_0 = Expression(('0.0', '20.0+273.15'), degree=1)
u_n = project(u_0, V)

g = Expression('10.0e3', degree=1)


# Boundary marker
class HeatFlux(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[0] > (l-DOLFIN_EPS))

# Initialize sub-domain instances
right = HeatFlux()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
right.mark(boundaries, 1)
dss = ds[boundaries]


# Define source term
f = Expression('0.0',degree=1)

# Define expressions used in variational forms
Dt  = Constant(dt)
k_A = Constant(k_A)
k2  = Constant(E_a/R)
n   = Constant(n)
rho = Constant(rho)
C_p = Constant(C_p)
k   = Constant(k)
H_r = Constant(H_r)


F = (alpha - alpha_1)*v_2*dx + (Dt*k_A*exp(-k2/u)*(1-alpha)**n)*v_2*dx + \
    (alpha - alpha_1)*v_2*dx - f*Dt*v_2*dx + \
    rho*C_p*u*v*dx + Dt*k*inner(nabla_grad(u), nabla_grad(v))*dx - (rho*C_p*u_1 - Dt*rho*H_r*f)*v*dx + Dt*g*v*dss(1)


T = 2.0 # Total simulation time: 8 minutes
t = 0.0
while t <= T:
    print 'time = %f' % t

    t += dt

    # Solve variational problem for time step
    solve(F == 0, u)

    # Update previous solution
    u_n.assign(u)


plot(u, title="Temperature")
interactive()

---

I tried adding a third trial function for the f = (alpha-alpha_1)/dt term and I got the same error