Function dependent on solution not updating each time step

variational formulation
1

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:

	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:

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.

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:


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))
2

changing the lines to

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 :slight_smile:

3

Thank you so mcuh! Still getting the hang of the ufl notation

4

Hello Emanuel
Is below lines alternative to your solution?

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
5

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 :slight_smile:

6

I will see that notation.
Thank you Emanuel