h.get_temp_map throwing value error

[hannan0786]

#Here are my material properties

prop_SiC = {} 
c_SiC = 10.08196 * u.angstrom
prop_SiC['heat_capacity'] = [24.298,
                              '63.38915 + (-1.23593e-4) * T - 0.39875e6 / T**2']
prop_SiC['therm_cond'] = [27.3,
                           1.48 * u.W / (u.m * u.K)]

g = 1.12e18

prop_SiC['sub_system_coupling'] = \
    ['-{:f} * (T_0 - T_1)'.format(g),
     '{:f} * (T_0 - T_1)'.format(g)]
prop_SiC['lin_therm_exp'] = [0, 4.5e-6]
prop_SiC['sound_vel'] = 8.3875 * u.nm / u.ps
prop_SiC['opt_ref_index'] = 2.553531 + 3.34j

SiC = ud.UnitCell('SiC', 'Silicon Carbide', c_SiC, **prop_SiC)
S = ud.Structure('Single Layer')
S.add_sub_structure(SiC, ### 100)

h = ud.HeatDiffusion()

then setting the following excitation parameters

h.excitation = {'fluence': [100] * u.mJ / u.cm**2,
                'delay_pump': [0] * u.ps,
                'pulse_width': [100] * u.ps,
                'multilayer_absorption': True,
                'wavelength': 343 * u.nm,
                'theta': 90 * u.deg}

use the following grid

delays = np.r_[-100:1000:0.005] * u.ps
init_temp = np.zeros_like(S.grid[0])

Calculating the temperature map is where am stuck

h.heat_diffusion = True
h.boundary_conditions = {'top_type': 'isolator', 'bottom_type': 'isolator'}
init_temp = np.ones([S.get_number_of_layers(), 2])
init_temp[:, 0] = 300
init_temp[:, 1] = 300
temp_map, delta_temp = h.get_temp_map(delays, init_temp)

Here am get the following value error: array must not contain infs or NaNs
Could you please assist me here. I am trying to move from COMSOl to this package as am extremely impressed by the chain of functionality it offers. It might be a silly mistake.

[dschick]

Hi @hannan0786 ,

Thanks for reporting this error and your willingness to move to open-source software :)

I think I could locate the origin of your error, and yes, it is kind of silly ...

In your example, you create a UnitCell for the Silicon Carbide, which is only necessary if you are going to do some hard X-ray scattering simulations at the end of your simulation chain assuming crystalline SiC.
If this is your actual goal, then you are missing one important step in building the UnitCell, namely adding some Atoms to it.
You can find an example here: https://udkm1dsim.readthedocs.io/en/latest/examples/structure.html#unit-cells

Without adding the Atoms the UnitCell has no mass, density, etc and hence the calculations fail.
Obviously, the silly thing about it is that you should get a related error warning instead of the current cryptic one.
See a new issue to solve that problem #125

Another approach is to use AmorphousLayers for which you explicitly have to provide the density for each layer.
You can do all transient simulations using this type of layer as well as dynamical magnetic scattering simulations XrayDynMag

Hope that solves your problem, but do not hesitate to ask for further help

Best

Daniel

[hannan0786]

Hi Daniel,
Thanks a lot for the help. I was expecting it to be silly error. I do need to simulate the behaviour in crystalline SiC and therefore need to create a unit cell. I had expected the mixed atom I create earlier was passed on to the structure.
I created the following unit cell

Si = ud.Atom('Si')
C = ud.Atom('C')
prop_SiC = {}
c_SiC = 10.08196*u.angstrom
prop_SiC['heat_capacity'] = [24.298,
                              '63.38915 + (-1.23593e-4)*T - 0.39875e6/T**2'
                             ]
prop_SiC['therm_cond'] = [27.3,
                           490*u.W/(u.m*u.K)]

g = 1.12e18

prop_SiC['sub_system_coupling'] = \
    ['-{:f}*(T_0-T_1)'.format(g),
     '{:f}*(T_0-T_1)'.format(g)
    ]
prop_SiC['lin_therm_exp'] = [0, 4.5e-6]
prop_SiC['sound_vel'] = 8.3875*u.nm/u.ps
prop_SiC['opt_ref_index'] = 2.553531 + 3.34j

SiC= ud.UnitCell('SiC', 'Silicon Carbide', c_SiC, **prop_SiC)
SiC.add_atom(Si, 0.5)
SiC.add_atom(C, 0.5)
S = ud.Structure('Single Layer')
S.add_sub_structure(SiC, 100)

I excite with the following laser pulse

h.excitation = {'fluence': [100]*u.mJ/u.cm**2,
                'delay_pump':  [0]*u.ps,
                'pulse_width':  [200]*u.fs,
                'multilayer_absorption': True,
                'wavelength': 1030*u.nm,
                'theta': 90*u.deg}
# temporal and spatial grid
delays = np.r_[-100:200:0.1]*u.ps
_, _, distances = S.get_distances_of_layers()

I am constantly been shown the following error

ValueError: cannot reshape array of size 255136 into shape (3002,112,2)
It seems to have a relation with my fluence. There is no error if i have very small values of fluences. But even in that case the results appear to be weird.

[dschick]

Hi @hannan0786 ,

I do need to simulate the behaviour in crystalline SiC and therefore need to create a unit cell.

Just to make things clear, the behaviour in the 2TM and even phonon dynamics does not differ if you define both the UnitCell and AmorpousLayer with the same thermo-elastic parameters. The only difference appears in the hard X-ray diffraction simulations XrayDyn. But anyways, since it makes no difference, you can keep it this way.

I had expected the mixed atom I create earlier was passed on to the structure.

Unfortunately, the code can not guess what you are aiming at ;)

ValueError: cannot reshape array of size 255136 into shape (3002,112,2)

I have reproduced your example and for the exact same parameters, I do not get the error message. However, when I increase the fluence to 1000 mJ/cm² I get the same message.
This is most likely an issue with the numerics of the heat diffusion solver when reaching too high temperatures.
I am still wondering why it is working on my machine with your parameters.
May I ask on which hardware you are working?

The results for your current parameters look okay to me, at least when you zoom into the delay range by a factor of 100:

But in general, the temperatures are very high, especially for the phonons this is not realistic.
I am not aware of many materials which survive 100 mJ/cm² photoexcitation.

Best

Daniel

[hannan0786]

Just to make things clear, the behaviour in the 2TM and even phonon dynamics does not differ if you define both the UnitCell and AmorpousLayer with the same thermo-elastic parameters

This is something which I didn't presume. For the two-temperature model in crystalline semiconductors should consider the crystal symmetry, anisotropy, and phonon dispersion relations. The last one I guessed being considered in how you have modelled it.

1000mJ/cm² is a very common laser fluence when working with high NA objectives. Far above the plastic deformation regime but an ordinary laser processing fluence.
Peak fluence

• maximal energy density per unit area of raw 4mm 200fs beam having 15µJ is 238.732 mJ/cm² and this before any focussing optics. The lens I am using is a High NA 20x micorscope objective. SO I don't understand the argument.

Unfortunately, the code can not guess what you are aiming at ;)

I messed up while reading how to create a mixed atom if it is not passed to the lattice where it is invoked later.

I still get the same error for anything above 30mJ/cm².
I haven't tried running this code on a cluster. I am using my Intel(R) Core(TM) i5-7500 CPU @ 3.40GHz with 16GB.

[dschick]

For the two-temperature model in crystalline semiconductors should consider the crystal symmetry, anisotropy, and phonon dispersion relations.

None of these properties is included explicitly. The simulation is 1D which excludes the first to properties.
The 2TM is completely built on top of the heat-diffusion equations on a continuous medium with quasi-equilibrium thermal properties.

• 4mm 200fs beam having 15µJ is 238.732 mJ/cm²

I think, here, we have a large discrepancy: the fluence in this code is simply defined as F = E/A where E is the pulse energy (in your case 15µJ) and A is the pulse area (pi (4mm/2)² (FWHM) in your case).
With these numbers, I end up with F = 0.119 mJ/cm² if I am not mistaken.
Could you please verify this. The pulse duration is not relevant here, but only for the 2TM.

However, I still do not get, why are receiving errors at 30mJ/cm² and I don't.
What are your python, numpy and scipy versions?
Could you send me your full code as a .py or .ipynb file?

[hannan0786]

Thanks Daniel. I was using a processing guy's Jargonised fluence term which is the peak fluence and is function of the Pulse Duration. As we want maximal energy density per unit area (at beam center).

I am starting to play around but it seem the errors now are more to deal with my system limitations if am not wrong.

```MemoryError:` Unable to allocate 916. MiB for an array with shape (30000000, 4) and data type float64

[dschick]

I am starting to play around but it seem the errors now are more to deal with my system limitations if am not wrong.

MemoryError: Unable to allocate 916. MiB for an array with shape (30000000, 4) and data type float64`

yes that seems like a problem with your available memory.
So either try to allocate and asign more memory to Python or reduce your number of delays and or numbers of layers

For the number of layers you can also reduce the spatial grid by the .distances property of the heat class.
See https://udkm1dsim.readthedocs.io/en/latest/examples/heat.html#heat-diffusion-parameters for details

[dschick]

I am closing this issue for now.

Do not hesitate to reopen it in case you still face any related problems.