Skip to content

Commit

Permalink
Merge pull request #5277 from fabidie/PRopenPMDforLASY
Browse files Browse the repository at this point in the history 
Enable complex input for FromOpenPMDPulse profile
  • Loading branch information
@ikbuibui
ikbuibui authored Feb 17, 2025
2 parents aed7105 + fb0e75d commit 17d9158
Show file tree
Hide file tree
Showing 4 changed files with 72 additions and 23 deletions.
4 changes: 3 additions & 1 deletion include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,10 @@ namespace picongpu
/** Datatype of the field record
*
* openPMD needs this information at compile time.
* @attention: supported are float, double and std::complex<double>.
* In case of complex input, the simulation uses only the real part.
*/
using dataType = float_64;
using dataType = double;

//! Defining the propagation ( = time) and polarisation axes of the input data
static constexpr char const* polarisationAxisOpenPMD = "x";
Expand Down
38 changes: 30 additions & 8 deletions include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
# include "picongpu/fields/incidentField/Functors.hpp"
# include "picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def"

# include <pmacc/math/Complex.hpp>
# include <pmacc/memory/buffers/HostDeviceBuffer.hpp>

# include <algorithm>
Expand Down Expand Up @@ -53,6 +54,7 @@
* - get rid of the 'wrong' transformation from time to space (z = c*t) of the
* longitudinal axis by using several iterations inside the openPMD file
* instead of just one
* - allow 2D simulations
*/

namespace picongpu
Expand All @@ -79,6 +81,7 @@ namespace picongpu
{
//! User SI parameters
using Params = T_Params;
using dataType = typename Params::dataType; // field record data type

//! Unit propagation direction vector in 3d
static constexpr float_X DIR_X = static_cast<float_X>(Params::DIRECTION_X);
Expand Down Expand Up @@ -149,8 +152,14 @@ namespace picongpu
or Params::polarisationAxisOpenPMD == "z"));

PMACC_CASSERT_MSG(
_error_propagationAxisOpenPMD_and_polarisationAxisOpenPMD_have_to_be_different_____check_your_parameters,
_error_propagationAxisOpenPMD_and_polarisationAxisOpenPMD_have_to_be_different____check_your_parameters,
(Params::polarisationAxisOpenPMD != Params::propagationAxisOpenPMD));

// check field data type
PMACC_CASSERT_MSG(
_error_unsupported_dataType____check_your_parameters,
(std::is_same<dataType, float>::value or std::is_same<dataType, double>::value
or std::is_same<dataType, std::complex<float>>::value));
};

template<typename T_Params>
Expand All @@ -172,6 +181,7 @@ namespace picongpu
{
//! Unitless parameters type
using Params = FromOpenPMDPulseUnitless<T_Params>;
//! field record data type
using dataType = typename Params::dataType;

//! FromOpenPMD pulse E functor
Expand Down Expand Up @@ -244,7 +254,7 @@ namespace picongpu
//! necessary attributes
// Raw = not yet aligned
::openPMD::Extent const extentRaw = meshRecord.getExtent();
auto const cellSizeRaw = mesh.gridSpacing<dataType>();
auto const cellSizeRaw = mesh.gridSpacing<float_X>();

bufferExtentOpenPMD = std::make_shared<pmacc::HostDeviceBuffer<float_X, 1u>>(3u);
bufferCellSizeOpenPMD = std::make_shared<pmacc::HostDeviceBuffer<float_X, 1u>>(3u);
Expand Down Expand Up @@ -304,18 +314,30 @@ namespace picongpu
openPMDIdx[aligningAxisIndex[d]] = tmpIndex % extentRaw[d];
tmpIndex /= extentRaw[d];
}
hostFieldDataBox(openPMDIdx)
= static_cast<float_X>(fieldData.get()[linearIdx] * meshRecord.unitSI())
/ sim.unit.eField();
if constexpr(std::is_same<dataType, std::complex<float>>::value)
{
// if the field record is complex, take only its real part
hostFieldDataBox(openPMDIdx)
= static_cast<float_X>(fieldData.get()[linearIdx].real() * meshRecord.unitSI())
/ sim.unit.eField();
}
else if constexpr(
std::is_same<dataType, float>::value or std::is_same<dataType, double>::value)
{
// field record is real
hostFieldDataBox(openPMDIdx)
= static_cast<float_X>(fieldData.get()[linearIdx] * meshRecord.unitSI())
/ sim.unit.eField();
}
}

/* If the transversal simulation window is smaller than the transversal DataBox extent,
* data at the chunk borders will be discarded. The maximum discarded value (relative to
* the maximum amplitude) will be logged.
*/
dataType maxE(0.0);
dataType maxEDiscarded(0.0);
dataType valE, valRight, valLeft;
float_X maxE(0.0_X);
float_X maxEDiscarded(0.0_X);
float_X valE, valRight, valLeft;
bool discard = false;
for(uint32_t d = 0u; d < 3u; d++)
{
Expand Down
45 changes: 34 additions & 11 deletions lib/python/picongpu/extra/input/preparingInsightData.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
import h5py
from scipy.interpolate import RegularGridInterpolator
from scipy.optimize import curve_fit
import scipy.constants as const

"""
ATTENTION!
Expand All @@ -29,7 +30,7 @@
"""

# please correct if other units than mm and fs are used!
c = 2.99792e-4 # speed of light in mm/fs
c = const.c * 1e-12 # speed of light in mm/fs


def gauss2D(xy, amp, x0, y0, w):
Expand Down Expand Up @@ -742,7 +743,9 @@ def to_time_domain(self, Ew, lamb_supp=10):

print("field data size: %.1f GB" % (np.prod(self.Et.shape) * 8 * 10**-9)) # assumes datatype double

def save_to_openPMD(self, outputpath, outputname, energy, pol="y", crop_x=0, crop_y=0, crop_t=0):
def save_to_openPMD(
self, outputpath, outputname, energy, pol="y", crop_x=0, crop_y=0, crop_t_neg=0, crop_t_pos=0, is_complex=False
):
"""
Save the field data in time domain to an openPMD file. This output will be used for the
FromOpenPMDPulse profile.
Expand All @@ -762,16 +765,38 @@ def save_to_openPMD(self, outputpath, outputname, energy, pol="y", crop_x=0, cro
energy: beam energy in Joule. Is used to determine the correct amplitude in the time domain.
For that, the approximation z = c*t is used, which holds only for tau << zR/c!
pol: polarisation direction, either "x" or "y". Default is "y" (the long axis).
crop_x/y/t: if the field data chunk is too big, one can crop the borders from both sides by the length
given with these values (same unit as the corresponding scales). Default is 0.
crop_x/y: if the field data chunk is too big, one can crop the spatial extent from both sides by this
length (in mm). Default is 0.
crop_t_pos/neg: if the field data chunk is too big, one can crop the temporal extent by this value
(in fs). Default is 0.
is_complex: whether Et should be stored as complex-valued (np.complex64). Default is False (i.e. only
the real part of Et is written to openPMD)
"""
assert pol == "x" or pol == "y", "Oops, invalid polarisation direction!"
idx_x = int(crop_x / self.dx + 0.5)
idx_y = int(crop_y / self.dy + 0.5)
idx_t = int(crop_t / self.dt + 0.5)
idx_t_neg = int(crop_t_neg / self.dt + 0.5)
idx_t_pos = int(crop_t_pos / self.dt + 0.5)
Nx, Ny, Nt = np.shape(self.Et)
# we only need the real part of the field
E_save = np.array(np.real(self.Et[idx_y : Ny - idx_y, idx_x : Nx - idx_x, idx_t : Nt - idx_t]))

if is_complex:
# store the complete complex field as np.complex64
# info: openPMD does not support np.complex128 yet
E_save = np.array(
self.Et[idx_y : Ny - idx_y, idx_x : Nx - idx_x, idx_t_neg : Nt - idx_t_pos], dtype=np.complex64
)
I_sum = np.sum(np.real(E_save) ** 2)
else:
# we only need the real part of the field
E_save = np.array(
np.real(self.Et[idx_y : Ny - idx_y, idx_x : Nx - idx_x, idx_t_neg : Nt - idx_t_pos]), dtype=np.float64
)
I_sum = np.sum(E_save**2)

# if the field was cropped, look how much pulse energy was discarded
if crop_x > 0 or crop_y > 0 or crop_t_neg > 0 or crop_t_pos > 0:
I_sum_orig = np.sum(np.real(self.Et))
print(f"Discarded pulse energy due to smaller window: {(1-I_sum/I_sum_orig):.2%}")

series = openpmd.Series(outputpath + outputname, openpmd.Access.create)
ite = series.iterations[0] # use the 0th iteration
Expand Down Expand Up @@ -813,7 +838,7 @@ def save_to_openPMD(self, outputpath, outputname, energy, pol="y", crop_x=0, cro
# B_z = -+ i/w * d/d_trans E_pol; but neglectable since B_z << B_trans
# 1st sign for pol = "x", 2nd for "y"
dV = self.dx * self.dy * self.dt * c * 10**-9 # m**3
W = dV * 8.854e-12 * np.sum(E_save**2)
W = dV * const.epsilon_0 * I_sum
# correct the value if working with an applied aperture
if self.is_masked:
energy *= self.field_ratio
Expand All @@ -834,6 +859,4 @@ def save_to_openPMD(self, outputpath, outputname, energy, pol="y", crop_x=0, cro

del series

print(
"data successfully saved, field data size: %.f MB" % (np.prod(E_save.shape) * 8 * 10**-6)
) # assumes datatype double
print("data successfully saved, field data size: %.f MB" % (np.prod(E_save.shape) * 8 * 10**-6))
8 changes: 5 additions & 3 deletions lib/python/picongpu/extra/input/preparingInsightData_example.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -490,7 +490,8 @@
"source": [
"## Save data to openPMD\n",
"The data is now nearly ready te be used as FromOpenPMDPulse input. The amplitude of the pulse in the time domain still has to be corrected (= scaled to the actual beam energy in Joule) before saving it to an openPMD file at the provided destination path.\n",
"Please pay attention to the size of the field data chunk: its real part will be stored on each used GPU as a whole, but their memory is limited. To reduce the chunk size, one can trim the edges by `crop_x`, `crop_y` (in mm) and `crop_t` (in fs)."
"Please pay attention to the size of the field data chunk: its real part will be stored on each used GPU as a whole, but their memory is limited. To reduce the chunk size, one can trim the edges by `crop_x`, `crop_y` (in mm), `crop_t_neg` and `crop_t_pos` (in fs).\n",
"One can choose to store the field data as real (default, i.e. `is_complex=False`) or complex (i.e. `is_complex=True`). The `FromOpenPMDPulse` profile can work with both, but will use only the real part of the field in case of complex input."
]
},
{
Expand All @@ -506,9 +507,10 @@
"pol = \"y\" # polarisation axis (can be 'x' or 'y')\n",
"crop_x = 0.3 # mm, trim transversal x axis\n",
"crop_y = 0.3 # mm, trim transversal y axis\n",
"crop_t = 100 # fs, trim time axis\n",
"crop_t_neg = 100 # fs, trim time axis from below\n",
"crop_t_pos = 100 # fs, trim time axis from above\n",
"\n",
"insight.save_to_openPMD(\"\", \"insightData_prepared%T.h5\", E, pol, crop_x, crop_y, crop_t)"
"insight.save_to_openPMD(\"\", \"insightData_prepared%T.bp5\", E, pol, crop_x, crop_y, crop_t_neg, crop_t_pos)"
]
},
{
Expand Down

0 comments on commit 17d9158

Please sign in to comment.