Implementation of linear interpolation within Fenics

Linear Algebra
1

Hi everyone,

I am currently trying to implement some kind of linear interpolation algorithm for a heat equation simulation. In this simulation, I want to assume that some constants vary with temperature and wish to interpolate this function for a large range of temperature values given some real physical data. Then, I would like to feed this interpolated function back into my solver so that the simulation can be solved accordingly.

For example if I assume I have some data here showing how a certain constant varies:

T   k
300   130
500   120
700   110
900   100
1100   90
1300   80

How would I interpolate this so that our constant, k , is defined over a large temperature range and then feed this back into my solver?
My heavily simplified example code is shown below:

from dolfin import *
from ufl import as_tensor

k = 130 #placeholder value for interpolation
# import data here????

mesh = UnitSquareMesh(10,10)
space = FunctionSpace(mesh, 'P', 1)
kappa_space = FunctionSpace(mesh, 'CG',1)
kappa = project(k, kappa_space)  #kappa is the tensor form of k

T = project(Expression("10*x[0]+100*x[1]", degree=2), Space)
print(kappa.vector().get_local())
2

What you would like to solve is a non-linear heat equation (as I assume k=kappa is the thermal conductivity), and the equation you are solving is \partial u/\partial t - \nabla \cdot (\kappa(u) \nabla u) = f.
Thus you should define \kappa as a function of T.
If this expression is not explicitly given, as in the case of the non-linear Poisson demo. You should create an analytical expression using for instance numpy’s polyfit function as illustrated below:

from dolfin import *
import numpy as np

k = [130,120,110,100,90,80]
T = [300,500,700,900,1100,1300]

weights = np.polyfit(T, k, deg=1)

kappa = lambda T: weights[0]*T+weights[1]

print(kappa(300), kappa(1500))

mesh  = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 1)
u,v = Function(V), TestFunction(V)
dt = 0.1
u_prev = Function(V)
x = SpatialCoordinate(mesh)
F = 1/dt*inner(u-u_prev, v)*dx + kappa(u)*inner(grad(u),grad(v))*dx \
    - x[0]*x[1]*v*dx
solve(F==0, u)
print(u.vector().get_local())