How to compute parameters in a lookup table using Fenics variables

General
1

Hi everyone,

I am trying to solve the two-dimensional heat conduction equation, considering the thermal conductivity as a function of temperature. I tried to use scipy.interpolate. to get the thermal conductivity for a specific temperature. However, I got an error because of the type of data. Any idea how I can solve this problem?

2

Without a minimal working example it is extremely difficult to provide help or advice. The solution for your problem likely requires a bespoke approach for your application.

3

Hi Nate,

The Python script is below. Note that this script was mainly generated by FEATool Multiphysics.

FEATool Multiphysics → FEniCS project conversion

File featool-fenics.py created 29-Apr-2025 11:15:36

try:
from fenics import *
except:
raise ImportError(“Could not find/import fenics.”)
from timeit import default_timer as timer
import atexit
comm = MPI.comm_world
size = comm.Get_size()
rank = comm.Get_rank()

import library

from scipy.interpolate import interp1d
import pandas as pd
import numpy as np

Functions

def tc_inter(T):
# read txt file
thermal_cond = pd.read_csv(“../Data/ThermalCond.txt”)

# temperature and thermal conductivity 
TV = thermal_cond.iloc[:,0].to_numpy()
TCv = thermal_cond.iloc[:,0].to_numpy()

# interpolate
f = interp1d(Tv, TCv,kind='cubic')

tc_val = f(T) 

return tc_val

Mesh and subdomains.

with HDF5File(comm, “/tmp/featool-fenics.h5”, “r”) as f:
mesh = Mesh()
f.read(mesh, “/mesh”, False)
subd = MeshFunction(“size_t”, mesh, mesh.topology().dim())
f.read(subd, “/mesh”)
bdr = MeshFunction(“size_t”, mesh, mesh.topology().dim()-1)
f.read(bdr, “/boundary”)
dx = Measure(“dx”, domain=mesh, subdomain_data=subd)
ds = Measure(“ds”, domain=mesh, subdomain_data=bdr)
n = FacetNormal(mesh)
nz = n[0]
nr = n[1]

Finite element spaces.

E1 = FiniteElement(“P”, mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, E1)

Function spaces.

T = Function(V)
T_t = TestFunction(V)
T_n = Function(V)

Model constants and expressions.

dt = 0.001
tmax = 0.2
t = 0.0
t_i = Expression(“t”, degree=1, t=0.0)
z = Expression(“x[0]”, degree=1)
r = Expression(“x[1]”, degree=1)
h_grid = CellDiameter(mesh) if “CellDiameter” in globals() else CellSize(mesh)

rho_ht_0 = Constant(0.5)
cp_ht_0 = Constant(7000)
u_ht_0 = Constant(0)
v_ht_0 = Constant(0)
k_ht_0 = tc_inter(T)
q_ht_0 = 200.01800.01800.0

Mass forms.

m = ((rrho_ht_0cp_ht_0)TT_t)*dx

Bilinear forms.

a = dt*(((rrho_ht_0cp_ht_0*u_ht_0)T.dx(0)T_t + (rk_ht_0)T.dx(0)T_t.dx(0) + (rrho_ht_0cp_ht_0v_ht_0)*T.dx(1)T_t + (rk_ht_0)*T.dx(1)*T_t.dx(1))*dx)

Linear forms.

f = ((dt*(rq_ht_0) + (rrho_ht_0*cp_ht_0)*T_n)*T_t)*dx

Boundary conditions.

dbc0 = DirichletBC(V, Constant(300), bdr, 1)
dbc1 = DirichletBC(V, Constant(300), bdr, 2)
dbc = [dbc0, dbc1]

Initial conditions.

dvar = [“T”]
T_n = project(Expression(“Tini(x[1])”, degree=1), V)
assign(T, T_n)

Solver parameters.

param = {
“newton_solver”: {
“linear_solver”: “default” if “default” in linear_solver_methods() else “default”,
“preconditioner”: “default” if “default” in krylov_solver_preconditioners() else “default”,
“maximum_iterations”: 20,
“relaxation_parameter”: 1.0,
“relative_tolerance”: 1e-06,
“absolute_tolerance”: 1e-06,
}
}

Solve.

f_data = HDF5File(comm, “/tmp/featool-fenics.h5”, “a”)
atexit.register(lambda: f_data.close())
t_start = timer()
f_data.write(T, “/” + dvar[0], 0.0)
n_ts = int(-(tmax // -dt))
for i_step in range(n_ts):
t += dt
print(“Time step %i, time = %f” % (i_step+1, t)) if rank == 0 else None
t_i.t = t
f = ((dt*(rq_ht_0) + (rrho_ht_0*cp_ht_0)*T_n)*T_t)*dx
solve(m + a - f == 0, T, dbc, solver_parameters=param)
T_n.assign(T)
f_data.write(T, “/” + dvar[0], t)
t_solve = comm.gather(timer() - t_start)
print(“t_solve = %g (tot: %g)” % (max(t_solve), sum(t_solve))) if rank == 0 else None