@VVKulagin the problem is that because of that division, if T=0 then the thermal conductivity is nan. It can be verified by adding model.log_level = 20.
One way to fix this is adding an epsilon to the division so that it’s never zero.
import festim as F
import numpy as np
import fenics as f
def thermal_cond_function(T):
return 1 + T**2 + 1 / (T + f.DOLFIN_EPS)
model = F.Simulation()
model.mesh = F.MeshFromVertices(vertices=np.linspace(0, 1, num=1000))
model.settings = F.Settings(
absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False
)
mat = F.Material(id=1, D_0=1, E_D=0.2, thermal_cond=thermal_cond_function)
model.materials = mat
model.T = F.HeatTransferProblem(transient=False, initial_condition=1)
model.boundary_conditions = [F.DirichletBC(surfaces=[1, 2], value=1 + F.x, field="T")]
model.initialise()
model.run()
Running the following code produces:
import matplotlib.pyplot as plt
f.plot(model.T.T)
plt.show()