From 226787433bb88bc7a83a3177495c2240237e857c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 23 Oct 2024 09:28:34 -0700 Subject: [PATCH 1/6] Add calculation of EM energy in RZ --- openpmd_viewer/addons/pic/lpa_diagnostics.py | 25 ++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/openpmd_viewer/addons/pic/lpa_diagnostics.py b/openpmd_viewer/addons/pic/lpa_diagnostics.py index 7f5d4d7..7b27b5a 100644 --- a/openpmd_viewer/addons/pic/lpa_diagnostics.py +++ b/openpmd_viewer/addons/pic/lpa_diagnostics.py @@ -675,6 +675,31 @@ def get_laser_envelope( self, t=None, iteration=None, pol=None, # Return the result return( envelope, info ) + + def get_electromagnetic_energy( self, t=None, iteration=None ): + """ + Compute the total electromagnetic energy inside the box (in Joules), i.e. + + .. math:: + + \\mathcal{E}_{tot} = \\int dV\,\\left[ \\frac{\\epsilon_0}{2}\\boldsymbol{E}^2 + \\frac{1}{2\mu_0}\\boldsymbol{B}^2 \\right] + + For simulations of high-intensity lasers propagating in underdense plasmas, + this is approximately equal to the laser energy (since the plasma fields are usually low) + """ + Ex, info = self.get_field('E', 'x', m=1, iteration=iteration ) + Ey, info = self.get_field('E', 'y', m=1, iteration=iteration ) + Ez, info = self.get_field('E', 'z', m=1, iteration=iteration ) + Bx, info = self.get_field('B', 'x', m=1, iteration=iteration ) + By, info = self.get_field('B', 'y', m=1, iteration=iteration ) + Bz, info = self.get_field('B', 'z', m=1, iteration=iteration ) + + energy_density = epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*mu_0)*(Bx**2 + By**2 + Bz**2) + volume = np.pi*abs(info.r)*info.dr*info.dz + E = (energy_density*volume[:, np.newaxis]).sum() + return E + + def get_main_frequency( self, t=None, iteration=None, pol=None, m='all', method='max'): """ From d312cdeb5a82a9f83d4f711af17b8260e2ae0ae5 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 23 Oct 2024 10:27:54 -0700 Subject: [PATCH 2/6] Update calculation --- openpmd_viewer/addons/pic/lpa_diagnostics.py | 35 +++++++++++++++----- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/openpmd_viewer/addons/pic/lpa_diagnostics.py b/openpmd_viewer/addons/pic/lpa_diagnostics.py index 7b27b5a..7aa996d 100644 --- a/openpmd_viewer/addons/pic/lpa_diagnostics.py +++ b/openpmd_viewer/addons/pic/lpa_diagnostics.py @@ -686,17 +686,34 @@ def get_electromagnetic_energy( self, t=None, iteration=None ): For simulations of high-intensity lasers propagating in underdense plasmas, this is approximately equal to the laser energy (since the plasma fields are usually low) - """ - Ex, info = self.get_field('E', 'x', m=1, iteration=iteration ) - Ey, info = self.get_field('E', 'y', m=1, iteration=iteration ) - Ez, info = self.get_field('E', 'z', m=1, iteration=iteration ) - Bx, info = self.get_field('B', 'x', m=1, iteration=iteration ) - By, info = self.get_field('B', 'y', m=1, iteration=iteration ) - Bz, info = self.get_field('B', 'z', m=1, iteration=iteration ) + Parameters: + ----------- + t : float (in seconds), optional + Time at which to obtain the data (if this does not correspond to + an existing iteration, the closest existing iteration will be used) + Either `t` or `iteration` should be given by the user. + + iteration : int + The iteration at which to obtain the data + Either `t` or `iteration` should be given by the user. + + Returns: + -------- + A float with the total electromagnetic energy inside the box (in Joules) + """ + # Get field array in full 3D + Ex, info = self.get_field('E', 'x', t=t, iteration=iteration, theta=None ) + Ey, info = self.get_field('E', 'y', m=1, iteration=iteration, theta=None ) + Ez, info = self.get_field('E', 'z', m=1, iteration=iteration, theta=None ) + Bx, info = self.get_field('B', 'x', m=1, iteration=iteration, theta=None ) + By, info = self.get_field('B', 'y', m=1, iteration=iteration, theta=None ) + Bz, info = self.get_field('B', 'z', m=1, iteration=iteration, theta=None ) + + # Compute the energy energy_density = epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*mu_0)*(Bx**2 + By**2 + Bz**2) - volume = np.pi*abs(info.r)*info.dr*info.dz - E = (energy_density*volume[:, np.newaxis]).sum() + volume = info.dx*info.dy*info.dz # Cell volume + E = (energy_density*volume).sum() return E From d0aa6b0cf05282cfcf7cfed2d32e68a4971b0e82 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 23 Oct 2024 10:54:16 -0700 Subject: [PATCH 3/6] Add automated test --- openpmd_viewer/addons/pic/lpa_diagnostics.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openpmd_viewer/addons/pic/lpa_diagnostics.py b/openpmd_viewer/addons/pic/lpa_diagnostics.py index 7aa996d..a884ee2 100644 --- a/openpmd_viewer/addons/pic/lpa_diagnostics.py +++ b/openpmd_viewer/addons/pic/lpa_diagnostics.py @@ -678,7 +678,7 @@ def get_laser_envelope( self, t=None, iteration=None, pol=None, def get_electromagnetic_energy( self, t=None, iteration=None ): """ - Compute the total electromagnetic energy inside the box (in Joules), i.e. + Compute the total electromagnetic energy inside the box, i.e. .. math:: @@ -702,7 +702,7 @@ def get_electromagnetic_energy( self, t=None, iteration=None ): -------- A float with the total electromagnetic energy inside the box (in Joules) """ - # Get field array in full 3D + # For RZ: use `theta=None` to get the full 3D field Ex, info = self.get_field('E', 'x', t=t, iteration=iteration, theta=None ) Ey, info = self.get_field('E', 'y', m=1, iteration=iteration, theta=None ) Ez, info = self.get_field('E', 'z', m=1, iteration=iteration, theta=None ) @@ -711,7 +711,7 @@ def get_electromagnetic_energy( self, t=None, iteration=None ): Bz, info = self.get_field('B', 'z', m=1, iteration=iteration, theta=None ) # Compute the energy - energy_density = epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*mu_0)*(Bx**2 + By**2 + Bz**2) + energy_density = const.epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*const.mu_0)*(Bx**2 + By**2 + Bz**2) volume = info.dx*info.dy*info.dz # Cell volume E = (energy_density*volume).sum() return E From 5771a8d19c33af19cd3d617528bdfbd932b952f7 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 23 Oct 2024 11:00:49 -0700 Subject: [PATCH 4/6] Add checks --- openpmd_viewer/addons/pic/lpa_diagnostics.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/openpmd_viewer/addons/pic/lpa_diagnostics.py b/openpmd_viewer/addons/pic/lpa_diagnostics.py index a884ee2..b1da837 100644 --- a/openpmd_viewer/addons/pic/lpa_diagnostics.py +++ b/openpmd_viewer/addons/pic/lpa_diagnostics.py @@ -702,13 +702,21 @@ def get_electromagnetic_energy( self, t=None, iteration=None ): -------- A float with the total electromagnetic energy inside the box (in Joules) """ + # Check that the required data is available + if 'E' not in self.avail_fields or 'B' not in self.avail_fields: + raise ValueError('The fields E and B are required to compute the electromagnetic energy.') + # Check that the fields have either the thetaMode or 3dcartesian geometry + if (self.fields_metadata['E']['geometry'] not in ['thetaMode', '3dcartesian']) or \ + (self.fields_metadata['B']['geometry'] not in ['thetaMode', '3dcartesian']): + raise ValueError('The electromagnetic energy can only be computed for 3D Cartesian and cylindrical simulations.') + # For RZ: use `theta=None` to get the full 3D field Ex, info = self.get_field('E', 'x', t=t, iteration=iteration, theta=None ) - Ey, info = self.get_field('E', 'y', m=1, iteration=iteration, theta=None ) - Ez, info = self.get_field('E', 'z', m=1, iteration=iteration, theta=None ) - Bx, info = self.get_field('B', 'x', m=1, iteration=iteration, theta=None ) - By, info = self.get_field('B', 'y', m=1, iteration=iteration, theta=None ) - Bz, info = self.get_field('B', 'z', m=1, iteration=iteration, theta=None ) + Ey, info = self.get_field('E', 'y', t=t, iteration=iteration, theta=None ) + Ez, info = self.get_field('E', 'z', t=t, iteration=iteration, theta=None ) + Bx, info = self.get_field('B', 'x', t=t, iteration=iteration, theta=None ) + By, info = self.get_field('B', 'y', t=t, iteration=iteration, theta=None ) + Bz, info = self.get_field('B', 'z', t=t, iteration=iteration, theta=None ) # Compute the energy energy_density = const.epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*const.mu_0)*(Bx**2 + By**2 + Bz**2) From 6129bb427c8721ec4b4cec85cd845f76e829c934 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 23 Oct 2024 11:06:31 -0700 Subject: [PATCH 5/6] Add automated test --- tests/test_laser_energy.py | 57 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 tests/test_laser_energy.py diff --git a/tests/test_laser_energy.py b/tests/test_laser_energy.py new file mode 100644 index 0000000..79f8c50 --- /dev/null +++ b/tests/test_laser_energy.py @@ -0,0 +1,57 @@ +""" +This test file is part of the openPMD-viewer. + +It checks that the function `get_electromagnetic_energy` works correctly + +Usage: +This file is meant to be run from the root directory of openPMD-viewer, +by any of the following commands +$ python tests/test_laser_energy.py +$ py.test +$ python -m pytest tests + +Copyright 2024, openPMD-viewer contributors +Authors: Remi Lehe +License: 3-Clause-BSD-LBNL +""" +from openpmd_viewer.addons import LpaDiagnostics +from scipy.constants import c +import numpy as np + +# Download required datasets +import os +def download_if_absent( dataset_name ): + "Function that downloads and decompress a chosen dataset" + if os.path.exists( dataset_name ) is False: + import wget, tarfile + tar_name = "%s.tar.gz" %dataset_name + url = "https://github.com/openPMD/openPMD-example-datasets/raw/draft/%s" %tar_name + wget.download(url, tar_name) + with tarfile.open( tar_name ) as tar_file: + tar_file.extractall() + os.remove( tar_name ) +download_if_absent( 'example-thetaMode' ) + +def test_laser_energy(): + """ + Check that the function `get_electromagnetic_energy` gives the same result + as the engineering formula for a Gaussian pulse + + E = 2.7e-5*(a0*w0/lambd)**2*(tau[fs]) + """ + ts = LpaDiagnostics('./example-thetaMode/hdf5') + iteration = 300 + + # Evaluate the laser energy using the engineering formula + a0 = ts.get_a0(iteration=iteration, pol='y') + w0 = ts.get_laser_waist(iteration=iteration, pol='y') + tau_fs = ts.get_ctau(iteration=iteration, pol='y') / c * 1e15 + omega = ts.get_main_frequency(iteration=iteration, pol='y') + lambd = 2*np.pi*c/omega + E_eng = 2.7e-5*(a0*w0/lambd)**2*tau_fs + + # Evaluate the laser energy using the function + E_func = ts.get_electromagnetic_energy(iteration=iteration) + + # Compare the two results + assert np.isclose(E_eng, E_func, rtol=0.04) From acb8d2566f83185513fb28e1f110aa43ad87a1d3 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 21 Feb 2025 14:21:59 -0800 Subject: [PATCH 6/6] Apply suggestions from code review --- openpmd_viewer/addons/pic/lpa_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openpmd_viewer/addons/pic/lpa_diagnostics.py b/openpmd_viewer/addons/pic/lpa_diagnostics.py index b1da837..4873fdb 100644 --- a/openpmd_viewer/addons/pic/lpa_diagnostics.py +++ b/openpmd_viewer/addons/pic/lpa_diagnostics.py @@ -721,7 +721,7 @@ def get_electromagnetic_energy( self, t=None, iteration=None ): # Compute the energy energy_density = const.epsilon_0/2.*(Ex**2 + Ey**2 + Ez**2) + 1./(2*const.mu_0)*(Bx**2 + By**2 + Bz**2) volume = info.dx*info.dy*info.dz # Cell volume - E = (energy_density*volume).sum() + E = energy_density.sum() * volume return E