Curve_fit of exponential function using FEM solution as input

I/O
1

Hi everyone,

I have fit an exponential function to some data, and now I want to use this exponential function with T, the solution to my finite element simulation, as an input. However, I am encountering the following error:

File "target-spotsize-fixedP.py", line 47, in exponenial_func
    return a*exp(-b/x)
ufl.log.UFLValueError: Invalid type conversion: [...] can not be converted to any UFL type

However I cannot see where my problem is, is there something wrong with using Project() in this case to output the value of this exponential function across my whole mesh?

Below is my simplified code:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from ufl import as_tensor
from ufl import Index
import math
from scipy.optimize import curve_fit
from mpi4py import MPI

mesh = UnitSquareMesh(16,16)

x = [0.01,10,20,30,40,50,60,80,100,1000,10000]  #T range
y = [0,0.1,0.36,0.61,0.72,0.78,0.82,0.84,0.87,0.99,0.99999] # y to fit


def exponenial_func(x, a, b):
    return a*np.exp(-b/x)
    
popt1, pcov1 = curve_fit(exponenial_func, x, y ,p0 = [2, 10])

Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x[0]+0.1*x[1]", degree=2), Space)


vp_space = FunctionSpace(mesh, 'CG',1)
vp_out = Function(vp_space, name = '$T$')
vp_ = exponenial_func(T,*popt1)
vp_out = project(vp_, vp_space,solver_type='cg')
2

exp is overloaded (where possible, in general never use from some_module import *).
Since curve_fit belongs to scipy, you want to use scipy.exp (or related, such as numpy.exp). At the moment your code is using ufl.exp. scipy doesn’t know about ufl.

3

Unfortunately using scipy.exp or np.exp gives me another error, namely:

 File "target-spotsize-fixedP.py", line 243, in <module>
    vp_ = exponenial_func(T,*popt1)
  File "target-spotsize-fixedP.py", line 48, in exponenial_func
    return a*scipy.exp(-b/x)
TypeError: loop of ufunc does not support argument 0 of type Division which has no callable exp method
4

You’ll need to provide a minimal working example of the problem. Your code does not work in any case,

NameError: name 'temp_vp' is not defined
5

Ah, thanks for pointing that out. I have altered the MWE accordingly.

6

It’s not really clear what you’re trying to do. You’re setting vp_ = exponenial_func(T,*popt1), but exponenial_func is a normal python function (taking “x” as first argument) and T is declared as an undefined dolfin Function, T = Function(Space). What could exponenial_func(T) possibly mean?

7

What I am trying to do here is take the solution T, from across my mesh and then calculate the value of exponenial_func according to the fit calculated earlier with the solution, T at a point on the mesh, taken as the x argument in the function. My hope would then this would calculate the value of exponenial_func across the whole mesh.

8

I think I see what you mean. But your code is not doing that. You are not providing any earlier solution to your exponenial_func. You’ve declared T as an empty (undefined) function, T = Function(Space) . It has no value.

Also your sample code does not run anyway (never gets to the line you’re trying to fix). Please make sure you’ve run and tested your code.

9

Ok, so I have updated the code slightly, now it is at the line that I wish to fix. I have also altered exp() to np.exp() but there is still this issue that I quoted in the original question.

10

Right, so now your dolfin function T is defined, and you want to pass it into your exponenial_func.

If I understand correctly, your situation is that you have some parametric function f(t) (your exponential function), and some other function T(x) (your dolfin function). What you’re trying to do is generate a third function in x,

g(x) = f( T(x) )

The title in your initial question was confusing, it sounded like you want to apply curve fitting to your calculated solution T(x).

Probably what you want to do is set up your transform function as an Expression, not a python function. Then you can project it like you did with your T function. But the Expression needs to take T as a parameter. That’s no longer simple with a simple in-line Expression, instead try a UserExpression

class exponenial_expr(UserExpression):
    def __init__(self, t, a, b, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.a = a
        self.b = b
    def eval(self, values, x):
        values[0] = exponenial_func( self.t(x), self.a, self.b )
    def value_shape(self):
        return ()

cf. Using functions in expressions

Then you can project,

vp = project( exponenial_expr( T,*popt1 , degree=2), Space)

or you might prefer to interpolate

vp_interp = interpolate( exponenial_expr( T,*popt1 , degree=2), Space)
1 Like
11

Wonderful! That works excellently, thank you very much!