'assemble cells' taking far too long on rectangular mesh?

Question

Hi there,

I'm having a small problem with the run times of my code. I'm trying to solve the equation:
d/dt(F(T) T) = nabla(kappa(T).nabla(T)) + S,
where F and kappa both depend on T and S is a constant. F and kappa are evaluated from data which is then interpolated.

Below is a minimal code, but i have not been able to give out the data values i have been using if that's ok.

I have run the 'list_timings' command and have found that 'Assemble cells' is taking about 20 seconds on average, repeated about 100 times per timestep???!! I'm not sure why this is as i am just using a standard rectangular mesh.... it is also taking about 10-11 newton iterations per timestep, whereas in matlab this used to run in under 4 iterations.

Is there something crucial I am doing wrong here? I'm a newbie to fenics so i'm not sure whether i've misunderstood what 'assemble cells' is actually doing.

thanks for your help in advance

from dolfin import *
import sys, numpy, time
import scipy.io
from scipy.interpolate import interp1d


def interpolate_table(Tvalue, data):
    ti = data[:,0] 
    fi = data[:,1]
    n = len(ti)
    retval=0

    for i in range(0,n-1):
        retval += conditional(And(ge(Tvalue,ti[i]), lt(Tvalue,ti[i+1])),
                              fi[i] + (fi[i+1]-fi[i]) / (ti[i+1]-ti[i]) * (Tvalue-ti[i]), 0)
    retval += conditional(le(Tvalue,ti[0]), fi[0], 0)
    retval += conditional(ge(Tvalue,ti[-1]), fi[-1], 0)

    return retval

#############################################################

xtop = 10000; xbottom = 0

mesh = RectangleMesh(xbottom, xbottom, xtop, xtop, 30, 10)

V = FunctionSpace(mesh, 'Lagrange', 1)

# Boundary conditions
tol = 1E-14
def top_boundary(x,on_boundary):
    return on_boundary and abs(x[0]-xtop)<tol
def bottom_boundary(x,on_boundary):
    return on_boundary and abs(x[0]-xbottom)<tol

BC_top = DirichletBC(V, Constant(10.0), top_boundary)
BC_bottom = DirichletBC(V, Constant(225.0), bottom_boundary)
bcs = [BC_bottom, BC_top]

T = Function(V)
v = TestFunction(V)
T_prev = Function(V)

class k_data(Expression):
    def eval(self, value, x):
        value[0] = interpolate_table(T(x), kappa_data)

class Fexp_data(Expression):
    def eval(self, value, x):
        value[0] = interpolate_table(T(x), F_data)

class Fprev_data(Expression):
    def eval(self, value, x):
        value[0] = interpolate_table(T_prev(x), F_data)

kappa = k_data()
F = Fexp_data()
F_prev = Fprev_data()

class IC(Expression):
    def eval(self,value,x):
                  value[0]=-215.0/10000.0*x[0]+225.0


T_prev = interpolate(IC(),V)

t_final = 1E11
dt = 1E10

# Assemble equations

S = Constant(0.8E-6)
time_part = ((F*T - F_prev*T_prev)/dt ) *v*dx
diff_part = inner(kappa*grad(T),grad(v))*dx
source_part = S*v*dx

# Stationary problem
Res = time_part + diff_part - source_part
J = derivative(Res,T)
problem = NonlinearVariationalProblem(Res, T, bcs, J)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters['newton_solver']
prm['absolute_tolerance']=1E-6
prm['relative_tolerance'] = 1E-6
prm['maximum_iterations'] = 100
prm['lu_solver']['report'] = True
prm['linear_solver'] = 'gmres'
prm['preconditioner'] = 'ilu'

t = dt
while t<=t_final:
    solver.solve()
    plot(T, rescale=False, title = 't=' + str(t))
    t += dt
    T_prev.assign(T)
    sol << (T,t)
    list_timings()

interactive()

umberto · Answer 1 · Nov 13, 2015

Hi jel2278,

My experience is that User-defined expressions by subclassing from the dolfin class Expression are extremely slow.

'assemble cells' taking far too long on rectangular mesh?

1 Answer