replace analytical integration with numerical

test/test_amorphous_layer.py

@@ -18,9 +18,7 @@ def test_amorphous_layer():
     assert al.name == 'Amorphous Layer'
     assert al.thickness == 2.86*u.angstrom
     assert al.heat_capacity[0](300) == 10
-    assert al.int_heat_capacity[0](300) == 3000
     assert al.lin_therm_exp[0](300) == 1e-6
-    assert al.int_lin_therm_exp[0](300) == 0.0003
     assert al.therm_cond[0](300) == 1
     assert al.opt_pen_depth == 11*u.nm
     assert al.sound_vel == 5*(u.nm/u.ps)
@@ -32,45 +30,30 @@ def test_amorphous_layer():
     assert al.heat_capacity_str == ['10', '1000']
     assert al.heat_capacity[0](300) == 10
     assert al.heat_capacity[1](300) == 1000
-    assert al.int_heat_capacity_str == ['10*T', '1000*T']
-    assert al.int_heat_capacity[0](300) == 3000
-    assert al.int_heat_capacity[1](300) == 300000
     assert al.therm_cond_str == ['10', '1000']
     assert al.therm_cond[0](300) == 10
     assert al.therm_cond[1](300) == 1000
     assert al.lin_therm_exp_str == ['10', '1000']
     assert al.lin_therm_exp[0](300) == 10
     assert al.lin_therm_exp[1](300) == 1000
-    assert al.int_lin_therm_exp_str == ['10*T', '1000*T']
-    assert al.int_lin_therm_exp[0](300) == 3000
-    assert al.int_lin_therm_exp[1](300) == 300000
     # test temperature-dependent parameters for str function input
     al.heat_capacity = ['10*T', 'exp(300-T)+300']
     al.therm_cond = ['10*T', 'exp(300-T)+300']
     al.lin_therm_exp = ['10*T', 'exp(300-T)+300']
     assert al.heat_capacity_str == ['10*T', 'exp(300-T)+300']
     assert al.heat_capacity[0](300) == 3000
     assert al.heat_capacity[1](300) == 301
-    assert al.int_heat_capacity_str == ['5*T**2', '300*T - exp(300 - T)']
-    assert al.int_heat_capacity[0](300) == 450000
-    assert al.int_heat_capacity[1](300) == 89999.0
     assert al.therm_cond_str == ['10*T', 'exp(300-T)+300']
     assert al.therm_cond[0](300) == 3000
     assert al.therm_cond[1](300) == 301
     assert al.lin_therm_exp_str == ['10*T', 'exp(300-T)+300']
     assert al.lin_therm_exp[0](300) == 3000
     assert al.lin_therm_exp[1](300) == 301
-    assert al.int_lin_therm_exp_str == ['5*T**2', '300*T - exp(300 - T)']
-    assert al.int_lin_therm_exp[0](300) == 450000
-    assert al.int_lin_therm_exp[1](300) == 89999.0
     # check backward compatibility
     al.heat_capacity = ['lambda T: 10*T', 'lambda T: exp(300-T)+300']
     assert al.heat_capacity_str == ['10*T', 'exp(300-T)+300']
     assert al.heat_capacity[0](300) == 3000
     assert al.heat_capacity[1](300) == 301
-    assert al.int_heat_capacity_str == ['5*T**2', '300*T - exp(300 - T)']
-    assert al.int_heat_capacity[0](300) == 450000
-    assert al.int_heat_capacity[1](300) == 89999.0
     # check subsystem temperatures
     al.therm_cond = ['10*T_0 + 30*T_1', 'exp(300-T_1)+300']
     assert al.therm_cond[0](np.array([300, 300])) == 12000

test/test_unitCell.py

@@ -22,9 +22,7 @@ def test_unit_cell():
     assert uc.b_axis == 2.86*u.angstrom
     assert uc.c_axis == 2.86*u.angstrom
     assert uc.heat_capacity[0](300) == 10
-    assert uc.int_heat_capacity[0](300) == 3000
     assert uc.lin_therm_exp[0](300) == 1e-6
-    assert uc.int_lin_therm_exp[0](300) == 0.0003
     assert uc.therm_cond[0](300) == 1
     assert uc.opt_pen_depth == 11*u.nm
     assert uc.sound_vel == 5*(u.nm/u.ps)
@@ -36,45 +34,30 @@ def test_unit_cell():
     assert uc.heat_capacity_str == ['10', '1000']
     assert uc.heat_capacity[0](300) == 10
     assert uc.heat_capacity[1](300) == 1000
-    assert uc.int_heat_capacity_str == ['10*T', '1000*T']
-    assert uc.int_heat_capacity[0](300) == 3000
-    assert uc.int_heat_capacity[1](300) == 300000
     assert uc.therm_cond_str == ['10', '1000']
     assert uc.therm_cond[0](300) == 10
     assert uc.therm_cond[1](300) == 1000
     assert uc.lin_therm_exp_str == ['10', '1000']
     assert uc.lin_therm_exp[0](300) == 10
     assert uc.lin_therm_exp[1](300) == 1000
-    assert uc.int_lin_therm_exp_str == ['10*T', '1000*T']
-    assert uc.int_lin_therm_exp[0](300) == 3000
-    assert uc.int_lin_therm_exp[1](300) == 300000
     # test temperature-dependent parameters for str function input
     uc.heat_capacity = ['10*T', 'exp(300-T)+300']
     uc.therm_cond = ['10*T', 'exp(300-T)+300']
     uc.lin_therm_exp = ['10*T', 'exp(300-T)+300']
     assert uc.heat_capacity_str == ['10*T', 'exp(300-T)+300']
     assert uc.heat_capacity[0](300) == 3000
     assert uc.heat_capacity[1](300) == 301
-    assert uc.int_heat_capacity_str == ['5*T**2', '300*T - exp(300 - T)']
-    assert uc.int_heat_capacity[0](300) == 450000
-    assert uc.int_heat_capacity[1](300) == 89999.0
     assert uc.therm_cond_str == ['10*T', 'exp(300-T)+300']
     assert uc.therm_cond[0](300) == 3000
     assert uc.therm_cond[1](300) == 301
     assert uc.lin_therm_exp_str == ['10*T', 'exp(300-T)+300']
     assert uc.lin_therm_exp[0](300) == 3000
     assert uc.lin_therm_exp[1](300) == 301
-    assert uc.int_lin_therm_exp_str == ['5*T**2', '300*T - exp(300 - T)']
-    assert uc.int_lin_therm_exp[0](300) == 450000
-    assert uc.int_lin_therm_exp[1](300) == 89999.0
     # check backward compatibility
     uc.heat_capacity = ['lambda T: 10*T', 'lambda T: exp(300-T)+300']
     assert uc.heat_capacity_str == ['10*T', 'exp(300-T)+300']
     assert uc.heat_capacity[0](300) == 3000
     assert uc.heat_capacity[1](300) == 301
-    assert uc.int_heat_capacity_str == ['5*T**2', '300*T - exp(300 - T)']
-    assert uc.int_heat_capacity[0](300) == 450000
-    assert uc.int_heat_capacity[1](300) == 89999.0
     # check subsystem temperatures
     uc.therm_cond = ['10*T_0 + 30*T_1', 'exp(300-T_1)+300']
     assert uc.therm_cond[0](np.array([300, 300])) == 12000

udkm1Dsim/simulations/heat.py

@@ -32,7 +32,7 @@
 import numpy as np
 from scipy.optimize import brentq
 from scipy.interpolate import interp2d
-from scipy.integrate import solve_ivp
+from scipy.integrate import solve_ivp, quad
 from time import time
 from os import path
 import warnings
@@ -677,7 +677,7 @@ def get_temperature_after_delta_excitation(self, fluence, init_temp, distances=[
         except AttributeError:
             pass
 
-        int_heat_capacities = self.S.get_layer_property_vector('_int_heat_capacity')
+        heat_capacities = self.S.get_layer_property_vector('_heat_capacity')
         thicknesses = self.S.get_layer_property_vector('_thickness')
         masses = self.S.get_layer_property_vector('_mass_unit_area')
         # masses are normalized to 1Ang^2
@@ -694,8 +694,9 @@ def get_temperature_after_delta_excitation(self, fluence, init_temp, distances=[
                 del_E = dalpha_dz[i]*E0*thicknesses[idx]
 
                 def fun(final_temp):
-                    return (masses[idx]*(int_heat_capacities[idx][0](final_temp)
-                                         - int_heat_capacities[idx][0](init_temp[i, 0]))
+                    return (masses[idx]*quad(heat_capacities[idx][0],
+                                             init_temp[i, 0],
+                                             final_temp)[0]
                             - del_E)
                 final_temp[i, 0] = brentq(fun, init_temp[i, 0], 1e5)
         delta_T = final_temp - init_temp  # this is the temperature change

udkm1Dsim/simulations/phonons.py

@@ -31,7 +31,7 @@
 import numpy as np
 from os import path
 from time import time
-from scipy.integrate import solve_ivp
+from scipy.integrate import solve_ivp, quad
 from tqdm.notebook import tqdm, trange
 
 
@@ -262,38 +262,30 @@ def calc_sticks_from_temp_map(self, temp_map, delta_temp_map):
         delta_temp_map = np.reshape(delta_temp_map, [M, L, K])
 
         thicknesses = self.S.get_layer_property_vector('_thickness')
-        int_lin_therm_exps = self.S.get_layer_property_vector('_int_lin_therm_exp')
+        lin_therm_exps = self.S.get_layer_property_vector('_lin_therm_exp')
 
-        # evaluated initial integrated linear thermal expansion from T1 to T2
-        int_alpha_T0 = np.zeros([L, K])
-        # evaluated integrated linear thermal expansion from T1 to T2
-        int_alpha_T = np.zeros([L, K])
         sticks = np.zeros([M, L])  # the sticks inserted in the unit cells
         sticks_sub_systems = np.zeros([M, L, K])  # the sticks for each thermodynamic subsystem
 
-        # calculate initial integrated linear thermal expansion from T1 to T2
-        # traverse subsystems
-        for ii in range(L):
-            for iii in range(K):
-                int_alpha_T0[ii, iii] = int_lin_therm_exps[ii][iii](temp_map[0, ii, iii]
-                                                                    - delta_temp_map[0, ii, iii])
-
         # calculate sticks for all subsystems for all delay steps
         # traverse time
         for i in range(M):
             if np.any(delta_temp_map[i, :]):  # there is a temperature change
-                # Calculate new sticks from the integrated linear thermal
+                # Calculate new sticks from integrating the linear thermal
                 # expansion from initial temperature to current temperature for
                 # each subsystem
                 # traverse subsystems
                 for ii in range(L):
                     for iii in range(K):
-                        int_alpha_T[ii, iii] = int_lin_therm_exps[ii][iii](temp_map[i, ii, iii])
+                        sticks_sub_systems[i, ii, iii] = \
+                            thicknesses[iii] * np.exp(quad(lin_therm_exps[ii][iii],
+                                                           temp_map[0, ii, iii] -
+                                                           delta_temp_map[0, ii, iii],
+                                                           temp_map[i, ii, iii])[0]) \
+                            - thicknesses[iii]
 
                 # calculate the length of the sticks of each subsystem and sum
                 # them up
-                sticks_sub_systems[i, :, :] = np.tile(thicknesses, (K, 1)).T \
-                    * np.exp(int_alpha_T-int_alpha_T0) - np.tile(thicknesses, (K, 1)).T
                 sticks[i, :] = np.sum(sticks_sub_systems[i, :, :], 1)
             else:  # no temperature change, so keep the current sticks
                 if i > 0:

udkm1Dsim/structures/layers.py

@@ -30,7 +30,7 @@
 from .. import u, Q_
 import numpy as np
 from inspect import isfunction
-from sympy import integrate, lambdify, symbols, symarray
+from sympy import lambdify, symbols, symarray
 from tabulate import tabulate
 
 
@@ -84,12 +84,8 @@ class Layer:
            conductivity [W/(m K)].
         lin_therm_exp (list[@lambda]): list of T-dependent linear thermal
            expansion coefficient (relative).
-        int_lin_therm_exp (list[@lambda]): list of T-dependent integrated
-           linear thermal expansion coefficient.
         heat_capacity (list[@lambda]): list of T-dependent heat capacity
            function [J/(kg K)].
-        int_heat_capacity (list[@lambda]): list of T-dependent integrated heat
-           capacity function.
         sub_system_coupling (list[@lambda]): list of of coupling functions of
            different subsystems [W/m³].
         num_sub_systems (int): number of subsystems for heat and phonons
@@ -176,6 +172,7 @@ def check_input(self, inputs):
         # traverse each list element and convert it to a function handle
         for input in inputs:
             T = symbols('T')
+            argument = T
             if isfunction(input):
                 raise ValueError('Please use string representation of function!')
             elif isinstance(input, str):
@@ -189,25 +186,23 @@ def check_input(self, inputs):
                     # check for presence of indexing and use symarray as argument
                     if '_' in input:
                         T = symarray('T', k)
-                        output.append(lambdify([T], input, modules='numpy'))
-                    else:
-                        output.append(lambdify(T, input, modules='numpy'))
-                    output_strs.append(input.strip())
+                        argument = [T]
                 except Exception as e:
                     print('String input for layer property ' + input + ' \
                         cannot be converted to function handle!')
                     print(e)
             elif isinstance(input, (int, float)):
-                output.append(lambdify(T, input, modules='numpy'))
-                output_strs.append(str(input))
+                # nothing to do here
+                pass
             elif isinstance(input, object):
-                output.append(lambdify(T, input.to_base_units().magnitude, modules='numpy'))
-                output_strs.append(str(input.to_base_units().magnitude))
+                input = input.to_base_units().magnitude
             else:
                 raise ValueError('Layer property input has to be a single or '
                                  'list of numerics, Quantities, or function handle strings '
                                  'which can be converted into a lambda function!')
 
+            output.append(lambdify(argument, input, modules=['numpy', 'scipy']))
+            output_strs.append(str(input).strip())
         return output, output_strs
 
     def get_property_dict(self, **kwargs):
@@ -228,10 +223,9 @@ def get_property_dict(self, **kwargs):
         properties_by_types = {'heat': ['_thickness', '_mass_unit_area', '_density',
                                         '_opt_pen_depth', 'opt_ref_index',
                                         'therm_cond_str', 'heat_capacity_str',
-                                        'int_heat_capacity_str', 'sub_system_coupling_str',
-                                        'num_sub_systems'],
-                               'phonon': ['num_sub_systems', 'int_lin_therm_exp_str', '_thickness',
-                                          '_mass_unit_area', 'spring_const', '_phonon_damping'],
+                                        'sub_system_coupling_str', 'num_sub_systems'],
+                               'phonon': ['num_sub_systems', '_thickness', '_mass_unit_area',
+                                          'spring_const', '_phonon_damping'],
                                'xray': ['num_atoms', '_area', '_mass', '_deb_wal_fac',
                                         '_thickness'],
                                'optical': ['_c_axis', '_opt_pen_depth', 'opt_ref_index',
@@ -405,10 +399,6 @@ def heat_capacity(self):
     def heat_capacity(self, heat_capacity):
         # (re)calculate the integrated heat capacity
         self._heat_capacity, self.heat_capacity_str = self.check_input(heat_capacity)
-        # delete last anti-derivative
-        self._int_heat_capacity = None
-        # recalculate the anti-derivative
-        self.int_heat_capacity
 
     @property
     def therm_cond(self):
@@ -418,34 +408,6 @@ def therm_cond(self):
     def therm_cond(self, therm_cond):
         self._therm_cond, self.therm_cond_str = self.check_input(therm_cond)
 
-    @property
-    def int_heat_capacity(self):
-        if hasattr(self, '_int_heat_capacity') and isinstance(self._int_heat_capacity, list):
-            return self._int_heat_capacity
-        else:
-            self._int_heat_capacity = []
-            self.int_heat_capacity_str = []
-            T = symbols('T')
-            try:
-                for hcs in self.heat_capacity_str:
-                    integral = integrate(hcs, T)
-                    self._int_heat_capacity.append(lambdify(T, integral, modules='numpy'))
-                    self.int_heat_capacity_str.append(str(integral))
-            except Exception as e:
-                print('The sympy integration did not work. You can set the '
-                      'analytical anti-derivative of the heat capacity '
-                      'of your layer as function str of the temperature '
-                      'T by typing layer.int_heat_capacity = \'c(T)\' '
-                      'where layer is the name of the layer object.')
-                print(e)
-
-        return self._int_heat_capacity
-
-    @int_heat_capacity.setter
-    def int_heat_capacity(self, int_heat_capacity):
-        self._int_heat_capacity, self.int_heat_capacity_str = self.check_input(
-                int_heat_capacity)
-
     @property
     def lin_therm_exp(self):
         return self._lin_therm_exp
@@ -454,38 +416,6 @@ def lin_therm_exp(self):
     def lin_therm_exp(self, lin_therm_exp):
         # (re)calculate the integrated linear thermal expansion coefficient
         self._lin_therm_exp, self.lin_therm_exp_str = self.check_input(lin_therm_exp)
-        # delete last anti-derivative
-        self._int_lin_therm_exp = None
-        # recalculate the anti-derivative
-        self.int_lin_therm_exp
-
-    @property
-    def int_lin_therm_exp(self):
-        if hasattr(self, '_int_lin_therm_exp') and isinstance(self._int_lin_therm_exp, list):
-            return self._int_lin_therm_exp
-        else:
-            self._int_lin_therm_exp = []
-            self.int_lin_therm_exp_str = []
-            T = symbols('T')
-            try:
-                for ltes in self.lin_therm_exp_str:
-                    integral = integrate(ltes, T)
-                    self._int_lin_therm_exp.append(lambdify(T, integral, modules='numpy'))
-                    self.int_lin_therm_exp_str.append(str(integral))
-            except Exception as e:
-                print('The sympy integration did not work. You can set the '
-                      'analytical anti-derivative of the linear thermal expansion '
-                      'of your unit cells as lambda function of the temperature '
-                      'T by typing layer.int_lin_therm_exp = \'c(T)\' '
-                      'where layer is the name of the layer object.')
-                print(e)
-
-        return self._int_lin_therm_exp
-
-    @int_lin_therm_exp.setter
-    def int_lin_therm_exp(self, int_lin_therm_exp):
-        self._int_lin_therm_exp, self.int_lin_therm_exp_str = self.check_input(
-                int_lin_therm_exp)
 
     @property
     def sub_system_coupling(self):
@@ -549,12 +479,8 @@ class AmorphousLayer(Layer):
            conductivity [W/(m K)].
         lin_therm_exp (list[@lambda]): list of T-dependent linear thermal
            expansion coefficient (relative).
-        int_lin_therm_exp (list[@lambda]): list of T-dependent integrated
-           linear thermal expansion coefficient.
         heat_capacity (list[@lambda]): list of T-dependent heat capacity
            function [J/(kg K)].
-        int_heat_capacity (list[@lambda]): list of T-dependent integrated heat
-           capacity function.
         sub_system_coupling (list[@lambda]): list of of coupling functions of
            different subsystems [W/m³].
         num_sub_systems (int): number of subsystems for heat and phonons
@@ -677,12 +603,8 @@ class UnitCell(Layer):
            conductivity [W/(m K)].
         lin_therm_exp (list[@lambda]): list of T-dependent linear thermal
            expansion coefficient (relative).
-        int_lin_therm_exp (list[@lambda]): list of T-dependent integrated
-           linear thermal expansion coefficient.
         heat_capacity (list[@lambda]): list of T-dependent heat capacity
            function [J/(kg K)].
-        int_heat_capacity (list[@lambda]): list of T-dependent integrated heat
-           capacity function.
         sub_system_coupling (list[@lambda]): list of of coupling functions of
            different subsystems [W/m³].
         num_sub_systems (int): number of subsystems for heat and phonons