# Using a temperature dependent function in the heat equation **URL:** \ **Category:** variational formulation\ **Created:** [February 22, 2023, 9:47am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489 "2023-02-22T09:47:04Z")\ **Posts on this page:** 6\ **Page:** 1
**Author:** ![CoE2845](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/coe2845/32/1210_2.png) [@CoE2845](https://fenicsproject.discourse.group/u/CoE2845)\ **Post date:** [February 22, 2023, 9:47am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489/1 "2023-02-22T09:47:04Z")
Hi everyone, I am currently trying to model the cooling of a cylinder via radiation as a function of time using the heat equation, where the emissivity of the material `epsilon_eff` is dependent on the temperature. However, `epsilon_eff` seems not to be updating for each time step. The basic code is below: ```auto from fenics import * from dolfin import * # Define the domain r1 = 1.0 # inner radius r2 = 2.0 # outer radius L = 3.0 # length N = 20 # number of mesh elements # Create the mesh mesh = CylinderMesh(Point(0, -L/2, 0), Point(0, L/2, 0), r1, r2, N) # Define the function space V = FunctionSpace(mesh, 'P', 1) # emissivity function emiss_in = [0.03,0.042,0.054,0.062,0.080,0.091,0.109,0.116,0.135,0.194,0.213,0.232,0.251] temp_emiss = [300,400,500,600,700,800,900,1000,1100,1600,1800,2000,2200] emiss_f = UnivariateSpline(temp_emiss, emiss_in, k = 2, ext = 3) # Define the initial temperature u0 = Expression('1000', degree=1) # Define the boundary conditions bc = DirichletBC(V, u0, 'on_boundary') # Define the thermal conductivity and radiation coefficient k = Constant(1.0) sigma = Constant(5.67e-8) # Define the variational problem u = TrialFunction(V) v = TestFunction(V) f = Constant(0.0) emiss_eff_array = emiss_f(u.vector().get_local()) emiss_space_eff = FunctionSpace(mesh, 'CG',1) epsilon_eff = Function(emiss_space_eff) epsilon_eff.vector()[:] = emiss_eff_array a = k * inner(grad(u), grad(v)) * dx L = f * v * dx - epsilon_eff * sigma * u**4 * v * ds # Define the time-stepping method T = 100.0 # total simulation time dt = 0.1 # time step t = 0.0 # initial time u = Function(V) # Create the file for storing the solution file = File('cooling.pvd') # Time-stepping loop while t < T: print(t) t += dt solve(a == L, u, bc) print(max(u.vector())) print(max(epsilon_eff.vector())) file << (u, t) ``` I think that the issue maybe with the line, `epsilon_eff.vector()[:] = emiss_eff_array` but I am unsure about how this can be fixed. Thanks in advance ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 22, 2023, 10:02am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489/2 "2023-02-22T10:02:16Z")
You would have to call ```python emiss_eff_array = emiss_f(u.vector().get_local()) epsilon_eff.vector()[:] = emiss_eff_array ``` inside the loop for it to update. DOLFIN cannot track changes to the vector (as in assigning data directly) and create a dependency tree. Thus you need to update this manually whenever `u` changes. ---
**Author:** ![CoE2845](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/coe2845/32/1210_2.png) [@CoE2845](https://fenicsproject.discourse.group/u/CoE2845)\ **Post date:** [February 22, 2023, 10:33am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489/3 "2023-02-22T10:33:21Z")
I see. When I copy said lines into the time loop, I get the following error: ```auto RuntimeError: *** ------------------------------------------------------------------------- *** DOLFIN encountered an error. If you are not able to resolve this issue *** using the information listed below, you can ask for help at *** *** fenics-support@googlegroups.com *** *** Remember to include the error message listed below and, if possible, *** include a *minimal* running example to reproduce the error. *** *** ------------------------------------------------------------------------- *** Error: Unable to solve linear system using PETSc Krylov solver. *** Reason: Solution failed to converge in 2 iterations (PETSc reason DIVERGED_INDEFINITE_PC, residual norm ||r|| = 2.786253e+00). *** Where: This error was encountered inside PETScKrylovSolver.cpp. *** Process: 0 *** *** DOLFIN version: 2019.1.0 *** Git changeset: ba376b6aebd7a9bc089be46b50bdd9f5c548fb91 *** ------------------------------------------------------------------------- ``` ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 22, 2023, 11:20am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489/4 "2023-02-22T11:20:51Z")
Then you should inspect the values you set to `epsilon_eff.vector()` (for instance write `epsilon_eff` to file, and check if the values are sane. ---
**Author:** ![CoE2845](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/coe2845/32/1210_2.png) [@CoE2845](https://fenicsproject.discourse.group/u/CoE2845)\ **Post date:** [February 22, 2023, 11:25am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489/5 "2023-02-22T11:25:39Z")
The values of `epsilon_eff` are sane. I placed the lines after `solve(a == L, u, bc)` so I wonder if that affects it. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 22, 2023, 11:43am UTC](https://fenicsproject.discourse.group/t/using-a-temperature-dependent-function-in-the-heat-equation/10489/6 "2023-02-22T11:43:02Z")
They should of course affect the solver. You wanted the `epsilon_eff` to update every time `u` changes, which is during solve. Thus after the solving the problem, you need to update `epsilon_eff`, such that this is reflected in the next solve.