replace analytical integration with numerical by dschick · Pull Request #110 · dschick/udkm1Dsim

dschick · 2022-06-10T13:34:42Z

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 · 2022-06-10T15:08:09Z

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

dschick · 2022-06-14T08:00:50Z

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 · 2022-06-14T11:29:28Z

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

dschick · 2022-06-14T11:45:22Z

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 · 2022-06-15T09:18:34Z

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 · 2022-06-15T10:01:59Z

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 · 2022-06-15T11:13:53Z

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 · 2022-06-15T11:21:37Z

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 · 2022-06-15T12:03:05Z

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 · 2022-06-15T16:33:27Z

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 · 2022-06-16T12:08:54Z

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 · 2022-06-16T13:13:49Z

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 · 2022-06-16T13:48:00Z

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.

codecov · 2026-07-20T10:15:50Z

Codecov Report

❌ Patch coverage is 75.00000% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 86.02%. Comparing base (e54254b) to head (28cc962).

Files with missing lines	Patch %	Lines
udkm1Dsim/simulations/heat.py	75.00%	1 Missing ⚠️
udkm1Dsim/simulations/phonons.py	66.66%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #110      +/-   ##
===========================================
- Coverage    86.13%   86.02%   -0.11%     
===========================================
  Files           10       10              
  Lines         3216     3170      -46     
===========================================
- Hits          2770     2727      -43     
+ Misses         446      443       -3

Files with missing lines	Coverage Δ
udkm1Dsim/structures/layers.py	`87.34% <100.00%> (-0.68%)`	⬇️
udkm1Dsim/simulations/heat.py	`79.92% <75.00%> (ø)`
udkm1Dsim/simulations/phonons.py	`86.59% <66.66%> (+0.06%)`	⬆️

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

codspeed-hq · 2026-07-20T10:16:42Z

Merging this PR will degrade performance by 84.38%

⚠️

Different runtime environments detected

Some benchmarks with significant performance changes were compared across different runtime environments,
which may affect the accuracy of the results.

Open the report in CodSpeed to investigate

❌ 1 regressed benchmark
✅ 12 untouched benchmarks
🆕 1 new benchmark
⏩ 6 skipped benchmarks¹

Warning

Please fix the performance issues or acknowledge them on CodSpeed.

Performance Changes

	Benchmark	`BASE`	`HEAD`	Efficiency
❌	`test_calc_energy_map`	95.6 ms	612.1 ms	-84.38%
🆕	`test_calc_sticks_from_temp_map`	N/A	5 s	N/A

Tip

Investigate this regression by commenting @codspeedbot fix this regression on this PR, or directly use the CodSpeed MCP with your agent.

_{Comparing quadintegration (28cc962) with develop (e54254b)}

6 benchmarks were skipped, so the baseline results were used instead. If they were deleted from the codebase, click here and archive them to remove them from the performance reports. ↩

dschick · 2026-07-20T13:52:20Z

I did some performance testing using numerical vs analytical integration for the heat_capacity and lin_therm_exp summarized in the table below:

heat_capacity:

Heat.get_temperature_after_delta_excitation() 4.7ms --> 4.6ms
Heat.calc_energy_map() 8.4ms --> 54.8ms

lin_therm_exp:
3. Phonon.calc_sticks_from_temp_map() 38.3ms --> 460.1ms

So for 1. it actually does not matter much, at least for constant heat capacities.
For 2. the increase is significant, but I guess bearable as one is not doing it too often
In case of 3. the change is way too large.

So in my view there are two option:

A. just keep it as it is and integrate functions by hand if necessary
B. only drop integration for heat_capacity and keep it for lin_therm_exp which is probably much less used as temperature-dependent quantity and if functions might be more monotonic.

Daniel Schick added 3 commits June 10, 2022 15:33

replace analytical integration with numerical

652ad0a

change analytical integration to numerical

c6bde59

remove analytical integrals from layer

11e668a

Daniel Schick added 2 commits June 10, 2022 17:09

fix flake8

664f463

fix unit tests

34f0482

Merge branch 'develop' into quadintegration

a95fe7f

Daniel Schick added 2 commits June 16, 2022 13:58

add scipy to lambdify modules argument

0ff9938

simplify check_input

03f87bd

dschick force-pushed the develop branch from 3b180d7 to 0a027e2 Compare June 13, 2023 22:24

dschick added 6 commits July 13, 2026 23:44

Merge branch 'develop' into quadintegration

2fe0608

change check_input back to develop state

cfe2d1a

remove int_heat_capacity & int_lin_expansion from tests

e3110ea

use quad in calc_energy_map

4a4a7cc

polishing

aee7e1a

Merge branch 'develop' into quadintegration

378c6b6

dschick added 4 commits July 20, 2026 13:39

Merge branch 'develop' into quadintegration

6cc1a5e

Merge branch 'develop' into quadintegration

e573492

fix becnhmark test_calc_sticks_from_temp_map

7e1a8bd

Merge branch 'develop' into quadintegration

28cc962

dschick mentioned this pull request Jul 22, 2026

in-place replacement of sympy's analytical integration with scipy.integrate.quad #198

Open

Conversation

dschick commented Jun 10, 2022

Uh oh!

dschick commented Jun 10, 2022

Uh oh!

dschick commented Jun 14, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dschick commented Jun 14, 2022

Uh oh!

dschick commented Jun 14, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dschick commented Jun 15, 2022

Uh oh!

mamattern commented Jun 15, 2022

Uh oh!

dschick commented Jun 15, 2022

Uh oh!

dschick commented Jun 15, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dschick commented Jun 15, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mamattern commented Jun 15, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dschick commented Jun 16, 2022

Uh oh!

mamattern commented Jun 16, 2022

Uh oh!

dschick commented Jun 16, 2022

Uh oh!

codecov Bot commented Jul 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

codspeed-hq Bot commented Jul 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Merging this PR will degrade performance by 84.38%

Performance Changes

Footnotes

Uh oh!

dschick commented Jul 20, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

dschick commented Jun 14, 2022 •

edited

Loading

dschick commented Jun 14, 2022 •

edited

Loading

dschick commented Jun 15, 2022 •

edited

Loading

dschick commented Jun 15, 2022 •

edited

Loading

mamattern commented Jun 15, 2022 •

edited

Loading

codecov Bot commented Jul 20, 2026 •

edited

Loading

codspeed-hq Bot commented Jul 20, 2026 •

edited

Loading