# Function dependent on solution not updating each time step **URL:** \ **Category:** variational formulation\ **Created:** [December 2, 2021, 1:08pm UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217 "2021-12-02T13:08:40Z")\ **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:** [December 2, 2021, 1:08pm UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217/1 "2021-12-02T13:08:40Z")
I currently have a function of my solution `T` across my 3D mesh that I want to update with every time step, such that as `T` changes and converges towards a steady state solution, this function also updates: ```auto def Evap(T): return np.exp(-1/T)*(1/T)**0.5 evap_new_array = Evap(T.vector().get_local()) evap_space = FunctionSpace(mesh, 'CG',1) evap_rate = Function(evap_space) evap_rate.vector()[:] = evap_new_array ``` The equation I am trying to simulate is a heat equation so `T` represents temperature here. To solve this time dependent problem, I have a time loop: ```auto t = 0 t_start = 0 t_end = 300 dtime = 1 for t in np.arange(t_start + dtime,t_end + dtime,dtime): print( "Time", t) solver.solve() assign(T0, T) print('Maximum Temperature:', max(T.vector())) print('Evaporated Energy: ', assemble(evap_rate*ds)) ``` Where `T0` is the solution at the previous time step. However if I run this, the ‘Evaporated Energy’ term does not seem to update, even though `T` is increasing. ```auto Time 1.0 Max Temperature: 949.7284928755515 Evaporated Energy: 0.3452573833397461 Time 1.3333333333333333 Max Temperature: 1149.6234203318184 Evaporated Energy: 0.3452573833397461 Time 1.6666666666666665 Max Temperature: 1349.5090677144087 Evaporated Energy: 0.3452573833397461 Time 1.9999999999999998 Max Temperature: 1549.3940380192917 Evaporated Energy: 0.3452573833397461 Time 2.3333333333333335 Max Temperature: 1749.2789564179793 Evaporated Energy: 0.3452573833397461 Time 2.6666666666666665 Max Temperature: 1949.1638705497824 Evaporated Energy: 0.3452573833397461 Time 3.0 Max Temperature: 2149.0487842980638 Evaporated Energy: 0.3452573833397461 Time 3.3333333333333335 Max Temperature: 2348.93369800837 Evaporated Energy: 0.3452573833397461 ``` I am unsure if I need to explicitly update `evap_rate` in the time loop or if there is something wrong with how I have constructed it. My simplified code is given below: ```auto from __future__ import print_function import numpy as np import scipy as scipy import matplotlib.pyplot as plt from dolfin import * import ufl from ufl import as_tensor from ufl import Index import math from scipy.optimize import curve_fit from mpi4py import MPI parameters["allow_extrapolation"] = True parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["optimize"] = True set_log_level(50) #---------------- t_start = 0.0 t_end = 20 nstep = 60 dtime = (t_end - t_start)/ nstep #time step #--------------------------------------- T_am = 300 #ambient vacuum temperature (K) T_a4 = T_am**4 #--------------------------------------- mesh = UnitCubeMesh(10,10,10) #------------ P_in = 100 I = P_in Space = FunctionSpace(mesh, 'P', 1) T = Function(Space) T0 = Function(Space) T_init = Expression('Tambient', degree=1, Tambient=300.) T = project(T_init, Space) assign(T0,T) #------------------------------------------------ # We should set up function spaces for each constant that will change as a function of temperature #------------------------------------------------ def Evap(T): #g/cm^2-sec return np.exp(-1/T)*(1/T)**0.5 evap_new_array = Evap(T.vector().get_local()) evap_space = FunctionSpace(mesh, 'CG',1) evap_rate = Function(evap_space) evap_rate.vector()[:] = evap_new_array #------------------------------------------- V = TestFunction(Space) # Test Function used for FEA dT = TrialFunction(Space) # Temperature Derivative i = Index() G = as_tensor(T.dx(i), (i)) #gradient of T G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step q = -G F = (1/dtime*(T-T0)*V - q[i]*V.dx(i)) * dx + (evap_rate)*V*ds - I*V*ds #final form to solve Gain = derivative(F, T, dT) # Gain will be used as the Jacobian required to determine the evolution of a dynamic system problem = NonlinearVariationalProblem(F, T, [], Gain) solver = NonlinearVariationalSolver(problem) solver.parameters["newton_solver"]["relative_tolerance"] = 1E-4 solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-3 solver.parameters["newton_solver"]["convergence_criterion"] = "residual" solver.parameters["newton_solver"]["error_on_nonconvergence"] = True solver.parameters["newton_solver"]["linear_solver"] = "cg" solver.parameters["newton_solver"]["maximum_iterations"] = 10000 solver.parameters["newton_solver"]["preconditioner"] = "hypre_euclid" solver.parameters["newton_solver"]["report"] = True solver.parameters["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = False solver.parameters["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1E-4 solver.parameters["newton_solver"]["krylov_solver"]["absolute_tolerance"] = 1E-3 solver.parameters["newton_solver"]["krylov_solver"]["monitor_convergence"] = False solver.parameters["newton_solver"]["krylov_solver"]["maximum_iterations"] = 100000 t = 0 for t in np.arange(t_start + dtime,t_end + dtime,dtime): print( "Time", t) solver.solve() assign(T0, T) print(max(T.vector())) print('Evaporated Energy: ', assemble(evap_rate*ds)) ``` ---
**Author:** ![Emanuel\_Gebauer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/emanuel_gebauer/32/3004_2.png) [@Emanuel\_Gebauer](https://fenicsproject.discourse.group/u/Emanuel_Gebauer)\ **Post date:** [December 2, 2021, 2:20pm UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217/2 "2021-12-02T14:20:27Z")
changing the lines to ```auto def Evap(T): #g/cm^2-sec return exp(-1/T)*(1/T)**0.5 evap_rate = Evap(T) #evap_space = FunctionSpace(mesh, 'CG',1) #evap_rate = Function(evap_space) #evap_rate.vector()[:] = evap_new_array ``` does the trick for me 🙂 ---
**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:** [December 2, 2021, 3:44pm UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217/3 "2021-12-02T15:44:39Z")
Thank you so mcuh! Still getting the hang of the `ufl` notation ---
**Author:** ![ilia145](https://avatars.discourse-cdn.com/v4/letter/i/ee7513/32.png) [@ilia145](https://fenicsproject.discourse.group/u/ilia145)\ **Post date:** [December 27, 2021, 7:29am UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217/4 "2021-12-27T07:29:34Z")
Hello Emanuel Is below lines alternative to your solution? ```auto def Evap(T): #g/cm^2-sec return np.exp(-1/T)*(1/T)**0.5 F = (1/dtime*(T-T0)*V - q[i]*V.dx(i)) * dx + Evap(T)*V*ds - I*V*ds #final form to solve ``` ---
**Author:** ![Emanuel\_Gebauer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/emanuel_gebauer/32/3004_2.png) [@Emanuel\_Gebauer](https://fenicsproject.discourse.group/u/Emanuel_Gebauer)\ **Post date:** [January 3, 2022, 8:28am UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217/5 "2022-01-03T08:28:21Z")
I think you could run into problems using the numpy package in `Evap`, as this is a nonlinear problem and fenics has to calculate the derivative of the equation (note that `T` is a function in the system). Try to stick to the ufl-notation: [Form language — Unified Form Language (UFL) 2021.1.0 documentation](https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html) 🙂 ---
**Author:** ![ilia145](https://avatars.discourse-cdn.com/v4/letter/i/ee7513/32.png) [@ilia145](https://fenicsproject.discourse.group/u/ilia145)\ **Post date:** [January 3, 2022, 11:28am UTC](https://fenicsproject.discourse.group/t/function-dependent-on-solution-not-updating-each-time-step/7217/6 "2022-01-03T11:28:57Z")
I will see that notation. Thank you Emanuel