Using numpy functions and updating with each time step

I/O
1

Hi everyone,

I currently have a function Evap(T) that is a function of my solution T, I wish to include this in my equation that is solved in a time dependently, and have the value of this function update for every time step. Currently I have this as the function I wish to use:

def vp_function(x): # x is the values of the solution at each node
	a = 5.73623
	b = 13204.109
	c = -24.306
	return 10**(a - (b/(x + c)))

tempspace = np.linspace(200,10000,10000)
vp_f = vp_function(tempspace)

mesh = UnitSquareMesh(10,10)
Space = FunctionSpace(mesh, 'P', 1) 
T = Function(Space)
vp_new_array = vp_f(T.vector().get_local())
vp_space = FunctionSpace(mesh, 'CG',1)
vp = Function(vp_space)
vp.vector()[:] = vp_new_array

def Evap(T):
	return vp*(1/T)

Where Evap(T) is the function I wish to use in my final form. However when I run what I have so farm I get this:

  File "test.py", line 54, in <module>
    vp_new_array = vp_f(T.vector().get_local())
TypeError: 'numpy.ndarray' object is not callable

Is there a way that I can set vp up such that it is usable in Fenics and its value will update for every time step, such that the value of Evap(T) will change due to its 1/T dependence and its implicit T dependence from vp ? My current code is 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
from scipy.interpolate import interp1d
from scipy import interp
from scipy.interpolate import UnivariateSpline

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
set_log_level(50)
hfont = {'fontname':'Liberation Sans'}

#----------------
t_start = 0.0
t_end = 100
nstep = 200
dtime = (t_end - t_start)/ nstep  #time step


#---------IMPORT EVAPORATION DATA--------------
def vp_function(x):
	a = 5.73623
	b = 13204.109
	c = -24.306
	return 10**(a - (b/(x + c)))

tempspace = np.linspace(200,10000,10000)
vp_f = vp_function(tempspace)


#---------------------------------------
mesh = UnitSquareMesh(10,10)
Space = FunctionSpace(mesh, 'P', 1) 
T = Function(Space)
T0 = Function(Space)
T_init = Expression('Tambient', degree=1, Tambient=1000.)
T = project(T_init, Space)
assign(T0,T)


vp_new_array = vp_f(T.vector().get_local())
vp_space = FunctionSpace(mesh, 'CG',1)
vp = Function(vp_space)
vp.vector()[:] = vp_new_array

def Evap(T):
	return vp*(1/T)

#-------------------------------------------

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(T)*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(assemble(Evap(T)*da))
	print(vp(max(T.vector())))
2

Did you mean to use

vp_new_array = vp_function(T.vector().get_local())

?

As the error states, vp_f is a numpy.ndarray, not a function.