From 652ad0aa99023e65781d13a8533427cb64f3e807 Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Fri, 10 Jun 2022 15:33:04 +0200 Subject: [PATCH 1/7] replace analytical integration with numerical --- udkm1Dsim/simulations/heat.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/udkm1Dsim/simulations/heat.py b/udkm1Dsim/simulations/heat.py index 1c0a5d92..130b6771 100644 --- a/udkm1Dsim/simulations/heat.py +++ b/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 From c6bde594a2c55e42d244c17336cf77189e6b602e Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Fri, 10 Jun 2022 16:00:33 +0200 Subject: [PATCH 2/7] change analytical integration to numerical --- udkm1Dsim/simulations/phonons.py | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/udkm1Dsim/simulations/phonons.py b/udkm1Dsim/simulations/phonons.py index 1e3db20d..64bf62c9 100644 --- a/udkm1Dsim/simulations/phonons.py +++ b/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,29 @@ 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: From 11e668a709330e586460d3dffe0d6d830c853bbc Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Fri, 10 Jun 2022 17:06:04 +0200 Subject: [PATCH 3/7] remove analytical integrals from layer --- udkm1Dsim/structures/layers.py | 83 ++-------------------------------- 1 file changed, 3 insertions(+), 80 deletions(-) diff --git a/udkm1Dsim/structures/layers.py b/udkm1Dsim/structures/layers.py index 26592575..c7742fe1 100644 --- a/udkm1Dsim/structures/layers.py +++ b/udkm1Dsim/structures/layers.py @@ -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 @@ -228,10 +224,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 +400,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 +409,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 +417,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 +480,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 +604,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 From 664f463c7fbc4ba359efa896f2a9b13a64bf1204 Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Fri, 10 Jun 2022 17:09:36 +0200 Subject: [PATCH 4/7] fix flake8 --- udkm1Dsim/simulations/phonons.py | 3 ++- udkm1Dsim/structures/layers.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/udkm1Dsim/simulations/phonons.py b/udkm1Dsim/simulations/phonons.py index 64bf62c9..2b430c74 100644 --- a/udkm1Dsim/simulations/phonons.py +++ b/udkm1Dsim/simulations/phonons.py @@ -279,7 +279,8 @@ def calc_sticks_from_temp_map(self, temp_map, delta_temp_map): for iii in range(K): 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[0, ii, iii] - + delta_temp_map[0, ii, iii], temp_map[i, ii, iii])[0]) \ - thicknesses[iii] diff --git a/udkm1Dsim/structures/layers.py b/udkm1Dsim/structures/layers.py index c7742fe1..890d9c41 100644 --- a/udkm1Dsim/structures/layers.py +++ b/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 From 34f0482d795087c44564d52c82bfa71e9e286ca9 Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Fri, 10 Jun 2022 17:11:22 +0200 Subject: [PATCH 5/7] fix unit tests --- test/test_amorphous_layer.py | 17 ----------------- test/test_unitCell.py | 17 ----------------- 2 files changed, 34 deletions(-) diff --git a/test/test_amorphous_layer.py b/test/test_amorphous_layer.py index a02126ac..a8fe7125 100644 --- a/test/test_amorphous_layer.py +++ b/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,18 +30,12 @@ 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'] @@ -51,26 +43,17 @@ def test_amorphous_layer(): 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 diff --git a/test/test_unitCell.py b/test/test_unitCell.py index e08b0640..2a578dc1 100644 --- a/test/test_unitCell.py +++ b/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,18 +34,12 @@ 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'] @@ -55,26 +47,17 @@ def test_unit_cell(): 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 From 0ff9938d9b034c589c44216b8eba9227b2771a77 Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Thu, 16 Jun 2022 13:58:12 +0200 Subject: [PATCH 6/7] add scipy to lambdify modules argument --- udkm1Dsim/structures/layers.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/udkm1Dsim/structures/layers.py b/udkm1Dsim/structures/layers.py index 890d9c41..7253a739 100644 --- a/udkm1Dsim/structures/layers.py +++ b/udkm1Dsim/structures/layers.py @@ -185,19 +185,20 @@ 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')) + output.append(lambdify([T], input, modules=['numpy', 'scipy'])) else: - output.append(lambdify(T, input, modules='numpy')) + output.append(lambdify(T, input, modules=['numpy', 'scipy'])) output_strs.append(input.strip()) 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.append(lambdify(T, input, modules=['numpy', 'scipy'])) output_strs.append(str(input)) elif isinstance(input, object): - output.append(lambdify(T, input.to_base_units().magnitude, modules='numpy')) + output.append(lambdify(T, input.to_base_units().magnitude, + modules=['numpy', 'scipy'])) output_strs.append(str(input.to_base_units().magnitude)) else: raise ValueError('Layer property input has to be a single or ' From 03f87bd06d6fa1f2aa4d9b50f15e42551db1c874 Mon Sep 17 00:00:00 2001 From: Daniel Schick Date: Thu, 16 Jun 2022 14:06:46 +0200 Subject: [PATCH 7/7] simplify check_input --- udkm1Dsim/structures/layers.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/udkm1Dsim/structures/layers.py b/udkm1Dsim/structures/layers.py index 7253a739..09947b29 100644 --- a/udkm1Dsim/structures/layers.py +++ b/udkm1Dsim/structures/layers.py @@ -172,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): @@ -185,26 +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', 'scipy'])) - else: - output.append(lambdify(T, input, modules=['numpy', 'scipy'])) - 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', 'scipy'])) - 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', 'scipy'])) - 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):