Problem in estimating energy

[p0six-00]

Hello everyone, I am recently using Smilei and there is one thing I am missing. I read in the documentation that for a 2D simulation environment the kinetic energy obtained from the Scalar module, as well as the total and electromagnetic energy, should be in K_r N_r L_r^D, where D = 2 in my case. So with SI units, the unit of measurement should be J/m. To see if I get results that are consistent with what I want to simulate, I have done several tests and I don't understand how I can convert this value J/m with the energy that I actually want to put into the system with the laser.
In some test I had simulated confirmed the intuitive idea I had, that is it was enough to just multiply that value by the size of the simulation window divided by the wavelength of the laser.

Subsequent tests, however, disproved this, as in the image below:

Where the yellow column represents the values I expect for laser energy, and the green column the maximum value measured after simulation with Smilei. The last column represents the ratio between the two quantities.

Here how I calculate the a0 value for simulating my laser:

from scipy.constants import *

wavelength = 1064 * 10**-9
nu = c / wavelength
omega = 2*pi*nu

waist = 1 * 10**-6 #m
impulse = 1 * 10**-15 #fs
En = 1*10**-3 #J

area = pi * waist
# NOTE: Here it should be waist**2, but one thing I noticed was that the dependence of simulated energy on waist was linear and not quadratic. Even if it were waist**2, I have the same problem, however.

I = En / impulse / area
E_0 = np.sqrt(2 * I / (epsilon_0 * c))

a0 = e * E_0 / (m_e * c * omega)

Here my namelist simulation code:

import math
from scipy.constants import *


wavelength = 1.064 #um
frequency = c / (wavelength * 10**-6) / 10**15#Hz


l_real = 2 * pi / wavelength #1 um
t_real = 2 * pi * frequency #1 fs

l0 = 2. * pi # laser wavelength [in code units]
t0 = 2. * pi                  # optical cycle
Lsim = [30*l0, 30*l0]  # length of the simulation
Tsim = 50.*t_real            # duration of the simulation
resx = 16.               # nb of cells in one laser wavelength
rest = 30 # nb of timesteps in one optical cycle 

Main(
    geometry = "2Dcartesian",
    
    interpolation_order = 2 ,
    
    cell_length = [l0/resx,l0/resx],
    grid_length  = Lsim,
    
    number_of_patches = [ 8, 8 ],
    
    timestep = t0/rest,
    simulation_time = Tsim,
     
    EM_boundary_conditions = [['PML']],
)

LaserGaussian2D(
    a0              = 1, # normalized amplitude
    omega           = 1.,
    focus           = [ 15* l_real, 15* l_real], # coordinates of laser focus
    waist           = 1.*l_real,
    incidence_angle = 0.,
    time_envelope   = ttrapezoidal(start=0., plateau = 5.*t_real)
)


DiagScalar(
    every = rest
)

DiagFields(
    every = rest,
    fields = ['Ex','Ey','Ez','Bx','By','Bz']
)

Here how I read the results with happi:

import happi
import math
from scipy.constants import *
import matplotlib.pyplot as plt
import numpy as np


laser_wavelength_um = 1.064 *10**-6 #m
omega_r_SI          = 2*pi*c/(laser_wavelength_um)

S = happi.Open("sim1/results/", reference_angular_frequency_SI=omega_r_SI)

S.Scalar("Uelm",  units=['um', 'fs', "mJ/um"]).plot()
plt.show(block=True)

Maybe I'm getting lost in a glass of water, but I just can't figure out if I'm actually simulating the laser under the conditions I want. As a proof, I was hoping to make up for it in the total energy fed into the system, but I can't figure out how I can bring the value estimated from the simulations into an energy value.

[mccoys]

Are you sure happi managed to do the conversion? Try to remove the units arguments and see if the value is the same

[p0six-00]

I think there is something wrong with it.
If I specify both angular frequency and units of measurement:

If I specify only units of measurement without angular frequency:

If I specify only angular frequency without units:

Code:

import happi
import math
from scipy.constants import *
import matplotlib.pyplot as plt
import numpy as np


laser_wavelength_um = 1.064 *10**-6 #m
omega_r_SI          = 2*pi*c/(laser_wavelength_um)

S = happi.Open("sim1/results/", reference_angular_frequency_SI=omega_r_SI)

S.Scalar("Uelm",  units=['um', 'fs', "mJ/um"]).plot()
plt.show(block=True)

Regardless, what would be the way to convert the simulated J/m energy into the actual input energy?

[mccoys]

You might have found a bug. I will try to test it soon.

The code does not provide mJ/um, but the units that you wrote above. You must express those in SI units to make the conversion

[p0six-00]

Okay thank you very much!

But the Scalar diagnostic should provide an energy per unit of length in the 2D case (as specified by https://smileipic.github.io/Smilei/Understand/units.html#quantities-integrated-over-the-grid).
Therefore, to obtain the value of energy fed into the system through the laser (in J), it should be sufficient to multiply that result with the size of the system itself, isn't it?

[mccoys]

Multiplying by the waist should be close enough. The general idea is that the simulation considers space uniform in z, but you know that it is not the case in reality. The simulation is close to reality only within a distance z of the order of the waist.

[mccoys]

So the trick, in your case, is that the quantity $K_r L_r^2 N_r$ does not depend on $\omega_r$, so the conversion can be done with or without providing a value for $\omega_r$.

Anyways, your laser is not well injected. Check fields diagnostics

[gannettim]

Thank you mmcoys for your patience, your reasoning is clear to me. I repeated a few simulations in a different machine and indeed the two values (wanted and simulated) differ by a factor proportional to twice the waist as you suggested.

For the last thing you said, I didn't understand. I used the example provided in the tutorial with a few minor changes.

[mccoys]

Plot the results of your field diagnostics. You will see that the laser is injected mostly before the start of the simulation