Using derivative and solve functions for function that is not a variational statment

Question

I am trying to write a solver that recreates the results for a polymer curing setup from here:
https://www.comsol.com/blogs/modeling-the-thermal-curing-process/

The system is two coupled time dependent differential equations that can be solved in the same time loop using an implicit scheme. I can solve the Poison equation like in the examples but I can't solve the one for the polymer curing rate.

I looked at the UFL manual and wrote the implicit scheme like this:

F = alpha - alpha_1 + dt*k_A*exp(-k2/u)*(1-alpha)**n
J = derivative(F, alpha)
solve(F == 0, alpha, J)

To use a Newton method. But when I use the derivative I get on the derivative line: AttributeError: 'Sum' object has no attribute 'arguments'

And without it I get on the solve line:
NotImplementedError: Wrong number or type of arguments for overloaded function 'la_solve'.

What do these errors mean? Is it because F is not a variational statement?

The full code is shown below:

from dolfin import *

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

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
T_ambient = Expression('20+273.15', degree=2)  #In K!!!


# Get mesh and define function spaces
mesh = IntervalMesh(numElems, 0.0,l)
V = FunctionSpace(mesh, 'Lagrange', 2)

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)

ds = Measure("ds")[boundaries]

#Insulation boundary given already and not needed

# Initial condition
u_1 = project(T_ambient, V)

g = Expression('Q_heat', degree=2, Q_heat=Constant(Q_heat))

dt = 1.0      # time step

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)


alpha_1 = project(Expression('0.0', degree=2), V)
alpha   = project(Expression('0.0', degree=2), V)

f = (alpha-alpha_1)/dt

a = rho*C_p*u*v*dx + dt*k*inner(nabla_grad(u), nabla_grad(v))*dx
L = (rho*C_p*u_1 - dt*rho*H_r*f)*v*dx + dt*g*v*ds(1)

A = assemble(a)   # assemble only once, before the time stepping
b = None          # necessary for memory saving assemeble call


# Compute solution
u = Function(V)   # the unknown at a new time level
T = 2.0           # total simulation time: 8 minutes
t = dt

k2 =  E_a/R


while t <= T:
    print 'time = %f' % t
    b = assemble(L, tensor=b)

    #bc.apply(A, b)
    solve(A, u.vector(), b)

    t += dt

    u_1.assign(u)

    F = alpha - alpha_1 + dt*k_A*exp(-k2/u)*(1-alpha)**n
    J = derivative(F, alpha)
    solve(F == 0, alpha, J)

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

    alpha_1.assign(alpha)

    #uArray = u.vector().array()


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

Answer 1

I realize I'm rather late, but I just had a similar problem, so here are my 2c:

The error is indeed due to F not being a form. You need to use the weak form of the equation, see e.g. here for an example with the heat equation. Eyeballing it, I'd say you want something like

F = inner(p-alpha_1, q)*dx - dt*inner(k_A*exp(-k2/u)*(1-p)**n, q)*dx
J = derivative(F, p)
solve(F == 0, alpha, J=J)

where

p = Function(V)
q = TestFunction(V)

Notice that p is a Function.

You should however take this with a huge pinch of salt: I just quickly tried it out and the NewtonSolver doesn't converge, even after reducing the timestep and increasing the number of iterations. This might be because F is just wrong or because of any number of other reasons, see e.g. this answer.

Comments

Yes you're right that is why it wasn't working. I got it working I can send you the code if you need it

---

Sure! If it's not too long you can paste it here so anyone coming to this question will see it. Glad it worked :)

---

Here is the code:

from dolfin import *

#Constants (all in SI units)
rho = 1200.0       #density
C_p = 1000.0       #specific heat
k   = 0.2          #thermal conductivity
H_r = 5.0e5        #total heat of reaction
k_A = 20.0e3       #frequency factor
E_a = 50.0e3       #activation energy
R   = 8.314        #universal gas constant
n   = 1.4          #order of reaction

l      = 5.0e-3    #length of 1D interval
nElems = 100       #number of finite elements
dt     = 1.0       #time step

Q_heat    = 10.0e3 #10 kW/m^2 for one m^2 which is 10 kW of heating power
T_ambient = 20.0+273.15  #In Kelvins


# Get mesh and define function spaces
mesh = IntervalMesh(nElems, 0.0,l)
V = FunctionSpace(mesh, 'Lagrange', 2)

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)
ds = ds(subdomain_data = boundaries)

# Initial conditions
u_1     = project(Expression('T_ambient', degree=2, T_ambient=Constant(T_ambient)), V)
alpha_1 = project(Expression('0.0', degree=2), V)

# Boundary condition
g = Expression('Q_heat', degree=2, Q_heat=Constant(Q_heat))


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
alpha = TrialFunction(V)
v_2   = TestFunction(V)

alpha = project(Expression('0.0', degree=2), V)

a = rho*C_p*u*v*dx + dt*k*inner(nabla_grad(u), nabla_grad(v))*dx
A = assemble(a)   # assemble only once, before the time stepping
b = None          # necessary for memory saving assemeble call


# Compute solution
u  = Function(V)  # the unknown at a new time level
T  = 60.0*1.8     # total simulation time: 8 minutes
t  = 0.0
k2 = E_a/R

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

    b = assemble(L, tensor=b)
    solve(A, u.vector(), b)
    u_1.assign(u)

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

    t += dt
    print 'time = %f' % t

plot(u,     title="Temperature", range_min=0.0, range_max=700.0)
plot(alpha, title="Alpha",       range_min=0.0, range_max=1.0  )
interactive()