Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Profile from openPMD file #151

Merged
merged 18 commits into from
Aug 1, 2023
Merged

Conversation

MaxThevenet
Copy link
Contributor

@MaxThevenet MaxThevenet commented Jul 2, 2023

This PR creates a new profile to read a laser pulse from an openPMD file. It derives from fromArrayProfile. Works both when the profile read is a laser envelope or the electric field of a laser pulse. in the latter case, the envelope is obtained with a Hilbert transform. Both were tested and give the same result as shown in #141.

Example

from lasy.profiles.from_openpmd_profile import FromOpenPMDProfile
from lasy.laser import Laser

# Create profile: one of
# - from FBPIC output
profile = FromOpenPMDProfile('lab_diags/hdf5/', 10, (0,1), 'E', 'y')
# - from LASY output
profile = FromOpenPMDProfile('./laserfbpic', 0, (1,0), 'laserEnvelope', envelope=True, prefix='laser3d')

# Create laser
dim = 'xyt'
lo = (-120e-6, -120e-6, -10e-15)
hi = (+120e-6, +120e-6, +60e-15)
npoints=(32,64,128)
laser = Laser(dim, lo, hi, npoints, profile)
laser.write_to_file('./laserfbpic/laser3d')

# Plot
from openpmd_viewer import OpenPMDTimeSeries
import numpy as np
import matplotlib.pyplot as plt
ts = OpenPMDTimeSeries('./laserfbpic')
F, m = ts.get_field(field='laserEnvelope')
F = np.abs(F)
plt.figure(figsize=(8,8))
plt.imshow(F[:,:,F.shape[2]//2], aspect='auto')
plt.colorbar()

To do (probably in a separate PR):

  • Clean up handling of polarisation.
  • Some functions could be interesting in general, and could be in tools or (better IMO) openPMD-viewer.
  • Make compliant with latest version of the standard.
  • (only if needed) support all axes ordering.
  • The phase retrieval is 1D, which may cause problems for more complex laser pulses, so this could be improved.

@MaxThevenet MaxThevenet added the new feature request or implementation of a new feature label Jul 2, 2023
@MaxThevenet MaxThevenet requested a review from RemiLehe July 5, 2023 11:58
@MaxThevenet
Copy link
Contributor Author

As suggested by @RemiLehe FromOpenPMDProfile now derives from FromArrayProfile.

lasy/profiles/from_openpmd_profile.py Outdated Show resolved Hide resolved

if "z" in m.axes.values():
# Flip to get complex envelope in t assuming z = -c*t
h = np.flip(h, axis=2)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should probably use converters

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A typical workflow:

  • Run a FBPIC simulation, get z output
  • Read in LASY, and do this dummy conversion t=-z/c (needed as lasy requires t representation)
  • Write envelope to file, t representation
  • Read in from another code, and if needed do the inverse conversion (z = - c*t). Alternatively, this and the step above could be simplified if LASY could output z representation directly, see More features to the z-t converters #155.

Points 2 and 4 are exactly symmetric, which is what we want, and that's easily done with this dummy convertor, so I think it should be like this. But the clean way would be to make this dummy convertor an option of the import_from_z convertor, that's for another PR (after #155 is addressed).


if "z" in m.axes.values():
dt = m.dz / c
t = (m.z - m.z[0]) / c
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here you implicitly set t[0]=0, while we usually center beam in time so t_peak=0

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also to be discussed in #155.

Comment on lines 103 to 104
np.diff(phase),
weights=0.5 * (np.abs(h[:, :, 1:]) + np.abs(h[:, :, :-1])),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here you could probably use np.gradient(phase, axis=-1) with weights=np.abs(h)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed

lasy/profiles/from_openpmd_profile.py Outdated Show resolved Hide resolved
Comment on lines 113 to 120
if wavelength is None:
s = io.Series(path + "/" + prefix + "_%T.h5", io.Access.read_only)
it = s.iterations[iteration]
omg0 = it.meshes["laserEnvelope"].get_attribute("angularFrequency")
wavelength_loc = 2 * np.pi * c / omg0
else:
wavelength_loc = wavelength
array = F
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe if wavelength is provided and differs from the one in the field, we should physically re-define the field with the new wavelength, e.g.

            array = F
            s = io.Series(path + "/" + prefix + "_%T.h5", io.Access.read_only)
            it = s.iterations[iteration]
            omg0 = it.meshes["laserEnvelope"].get_attribute("angularFrequency")

            if wavelength is None:
                wavelength_loc = 2 * np.pi * c / omg0
            else:
                wavelength_loc = wavelength
                array *= np.exp(1j * (2 * np.pi * c  / wavelength - omg0) * t)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, we could do sth like this, I would keep this for another PR, as there are quite a few options that could make sense. I put it int #156.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i'm not sure it's an separate option -- imagine, a user takes an envelope defined with some central frequency, but wants to obtain an object with a different one. In that case the envelope phase should be adjusted as suggested. Otherwise wavelength_loc = wavelength is only passed to super().__init__() which registers it but never uses.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

likewise, discussed offline, this option is useful and independent from this PR.

@MaxThevenet MaxThevenet mentioned this pull request Jul 31, 2023
3 tasks
@MaxThevenet
Copy link
Contributor Author

Thanks for your comments! I should have addressed them. I suggest that, for some of them, we address them in a separate PR, so I opened an issue (#156) accordingly. In particular, there should be CI, but that would be vastly simplified by either openPMD supporting t representation or by working a bit on the converter (also related to #155).

@MaxThevenet
Copy link
Contributor Author

Discussed offline: I'll remove the wavelength argument, as it is mostly here for debugging and for some corner cases, and implement some other utils (find central frequency, shift the frequency at which the envelope is defined, etc.) that we can call here.

@MaxThevenet MaxThevenet requested review from hightower8083 and removed request for RemiLehe August 1, 2023 13:55
@MaxThevenet MaxThevenet merged commit b744334 into LASY-org:development Aug 1, 2023
4 checks passed
@MaxThevenet MaxThevenet deleted the from_openpmd branch August 1, 2023 14:00
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
new feature request or implementation of a new feature
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants