replace analytical integration with numerical

[dschick]

following #109 the analytical integration can be replaced by numerical integration using scipy.integrate.quad.
This avoids doing the symbolig integration of the heat capacities and lin_therm_exp coefficients

[dschick]

@mamattern, could you please check this branch, if it still works without analytical integration?

[dschick]

Here are some comparisons for different types of temperature dependence of the heat capacity using scipy.integrate.quad's default setting for the numerical integration.
The integrator keeps the rel. error well below the default value of 1.49e-8.
This should be fine for most of the use cases.
The rel. error of the integrator is compared to the rel. difference between the analytical and numerical integrations

[dschick]

Here I plot the same functions sampled in dT = 20K steps using np.interp

[dschick]

Following this thread quad encounters some problems if the function has any edges, which can be the result of a linear interpolation, e.g. from np.interp.
The proposed solution in the thread renders the integration rather slow.

Using scipy.interpolation.interp1d with quadratic interpolation results in a much better results, although the input data again sampled with dT=20K

In general the interpolation of experimental data should be avoided and a fit function should be used instead.
Since one would not need to analytically integrate the fit function of the heat_capacity and lin_therm_exp these could be more complex functions.

[dschick]

instead of using np.interp it is possible to use step-wise defined functions such as

def my_heat_capacity(T):
    T = np.array(T)
    res = np.ones_like(T, dtype=float)
    select = T>=634
    res[select] = -0.775/0.0718*np.abs((T[select]-634)/634)**0.0718 + 13.5+3.03*(T[select]-634)/634
    select = T<634
    res[select] =-2.46/0.0718*np.abs((T[select]-634)/634)**0.0718+37+3.03*(T[select]-634)/634
    return res

which results in a reasonable numerical integration

[mamattern]

Using the numerical integration branch and the above defined heat capacity by:

def my_heat_capacity(T):
    T = np.array(T)
    res = np.ones_like(T, dtype=float)
    select = T >= 634
    res[select] = (-0.775/0.0718*np.abs((T[select]-634)/634)**0.0718 + 13.5+3.03*(T[select]-634)/634)/58.69e-3
    select = T < 634
    res[select] = (-2.46/0.0718*np.abs((T[select]-634)/634)**0.0718+37+3.03*(T[select]-634)/634)/58.69e-3
    return res

prop_unit_cell['heat_capacity'] = ['lambda T: my_heat_capacity(T)', 400]

does not work and produces the following error:

  File "udkm1dsim_ni20nm_2tm_new.py", line 207, in <module>
    temp_map, delta_temp_map = h.get_temp_map(delays, Init_temp)

  File "ib\site-packages\udkm1Dsim\simulations\heat.py", line 740, in get_temp_map
    self.calc_temp_map(delays, init_temp)

  File "lib\site-packages\udkm1Dsim\simulations\heat.py", line 871, in calc_temp_map
    temp, _ = self.get_temperature_after_delta_excitation(fluence,

  File "lib\site-packages\udkm1Dsim\simulations\heat.py", line 701, in get_temperature_after_delta_excitation
    final_temp[i, 0] = brentq(fun, init_temp[i, 0], 1e5)

  File "lib\site-packages\scipy\optimize\_zeros_py.py", line 783, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)

  File "lib\site-packages\udkm1Dsim\simulations\heat.py", line 697, in fun
    return (masses[idx]*quad(heat_capacities[idx][0],

  File "lib\site-packages\scipy\integrate\_quadpack_py.py", line 351, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,

  File "lib\site-packages\scipy\integrate\_quadpack_py.py", line 463, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

  File "<lambdifygenerated-25>", line 2, in _lambdifygenerated
    [T_0, T_1] = _Dummy_26

[dschick]

the problem is raised by sympy's lambdify function which cannot handle the function my_heat_capacity
the quick workaround is to set the lambda function again directly via the private property:
Ni_layer._heat_capacity = [lambda T: my_heat_capacity(T)]
This should work, however, I noticed, that the calculation of the get_temperature_after_delta_excitation() took like 30sec for this complex heat_capacity. The same slow-down should also happen for the heat_diffusion in case the ode-solver is working with the heat_capacity

[dschick]

removing the vectorization of the heat_capacity function brings a huge speed improvement

def my_heat_capacity(t):
    if t >= 634:
        return -0.775/0.0718*np.abs((t-634)/634)**0.0718 + 13.5+3.03*(t-634)/634
    else:
        return -2.46/0.0718*np.abs((t-634)/634)**0.0718+37+3.03*(t-634)/634

[dschick]

another option with vectorization would be using scipy.special.erf function:

Ni_layer.heat_capacity = '(0.5*erf(634-T)+1)*(-2.46/0.0718*abs((T-634)/634)**0.0718+37+3.03*(T-634)/634) + (0.5*erf(T-634)+1)*(-0.775/0.0718*abs((T-634)/634)**0.0718 + 13.5+3.03*(T-634)/634)'

this again works with the lambdification and seems to be reasonable fast.

[mamattern]

This does not work for me. I import from scipy.special import erf and use the line above. However, python return an error:

 File "lib\site-packages\scipy\integrate\_quadpack_py.py", line 463, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

  File "<lambdifygenerated-1559>", line 2, in _lambdifygenerated
    return  (0.5*(erf(634-T)+1))*(-2.46/0.0718*abs((T-634)/634)**0.0718+37+3.03*(T-634)/634)/58.69e-3 + (0.5*(erf(T-634)+1))*(-0.775/0.0718*abs((T-634)/634)**0.0718 + 13.5+3.03*(T-634)/634)/58.69e-3

NameError: name 'erf' is not defined

If I use the private setting: Ni_layer._heat_capacity = [lamda T: ... erf(634-T) ..., 442.1] everything works fast and fine using the numeric integration branch.

[dschick]

I just fixed the errors from the lambdification.
@mamattern could you please try again?
Could you compare the current version with the original np.interp version in terms of accuracy and speed?

Would be very interesting for the phonon calculation as well

[mamattern]

Now it works for me. I tested for numeric integration both possibilities:

Ni_layer.heat_capacity = '(0.5*erf(634-T)+1)*(-2.46/0.0718*abs((T-634)/634)**0.0718+37+3.03*(T-634)/634) + (0.5*erf(T-634)+1)*(-0.775/0.0718*abs((T-634)/634)**0.0718 + 13.5+3.03*(T-634)/634)'

and

Ni_layer._heat_capacity = [lamda T: ... erf(634-T) ..., 442.1]

the difference in the results is of the order e-10 and the calculation of the temperature after delta excitation took 7.4/7.2s, the calculation of the heat diffusion took 189/179s and the calculation of the strain map took 3.0/2.9s.

As a comparison I also tested the analytic integration possibility by giving the heat capacity and linear thermal expansion coefficient and their integrals by np.interp(). This took 0.02s, 129s and 0.4s. However, the difference in the order of 5-10% in the strain if the divergence of the quantities is relevant. For temperatures well below, the difference is of the order of e-4. In contrast, the temperature (both phonon and electron) shows a difference of e-4 in both cases, so is significantly less affected.

[dschick]

okay, so without referring to the large difference in the strain_map (which I would still need to evaluate), I realize that the decrease in performance is not really balancing the increase in usability and code cleaning.

To that I end, I would stick to the analytical integration of the heat_capacity and lin_therm_exp for now. But I am thinking of maybe offering also the numerical integration as an option.
In addition, the option to manually set the analytical integral (e.g. by np.interp) should be better documented and maybe also offered in the according warning messages.