diff --git a/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def b/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def index a0e90f3f7f..86a170ce2f 100644 --- a/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def +++ b/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def @@ -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. + * 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"; diff --git a/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.hpp b/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.hpp index c9e04c308a..6cbb610bfe 100644 --- a/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.hpp +++ b/include/picongpu/fields/incidentField/profiles/FromOpenPMDPulse.hpp @@ -25,6 +25,7 @@ # include "picongpu/fields/incidentField/Functors.hpp" # include "picongpu/fields/incidentField/profiles/FromOpenPMDPulse.def" +# include # include # include @@ -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 @@ -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(Params::DIRECTION_X); @@ -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::value or std::is_same::value + or std::is_same>::value)); }; template @@ -172,6 +181,7 @@ namespace picongpu { //! Unitless parameters type using Params = FromOpenPMDPulseUnitless; + //! field record data type using dataType = typename Params::dataType; //! FromOpenPMD pulse E functor @@ -244,7 +254,7 @@ namespace picongpu //! necessary attributes // Raw = not yet aligned ::openPMD::Extent const extentRaw = meshRecord.getExtent(); - auto const cellSizeRaw = mesh.gridSpacing(); + auto const cellSizeRaw = mesh.gridSpacing(); bufferExtentOpenPMD = std::make_shared>(3u); bufferCellSizeOpenPMD = std::make_shared>(3u); @@ -304,18 +314,30 @@ namespace picongpu openPMDIdx[aligningAxisIndex[d]] = tmpIndex % extentRaw[d]; tmpIndex /= extentRaw[d]; } - hostFieldDataBox(openPMDIdx) - = static_cast(fieldData.get()[linearIdx] * meshRecord.unitSI()) - / sim.unit.eField(); + if constexpr(std::is_same>::value) + { + // if the field record is complex, take only its real part + hostFieldDataBox(openPMDIdx) + = static_cast(fieldData.get()[linearIdx].real() * meshRecord.unitSI()) + / sim.unit.eField(); + } + else if constexpr( + std::is_same::value or std::is_same::value) + { + // field record is real + hostFieldDataBox(openPMDIdx) + = static_cast(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++) { diff --git a/lib/python/picongpu/extra/input/preparingInsightData.py b/lib/python/picongpu/extra/input/preparingInsightData.py index ffd8790c6a..cd11f02d1f 100644 --- a/lib/python/picongpu/extra/input/preparingInsightData.py +++ b/lib/python/picongpu/extra/input/preparingInsightData.py @@ -9,6 +9,7 @@ import h5py from scipy.interpolate import RegularGridInterpolator from scipy.optimize import curve_fit +import scipy.constants as const """ ATTENTION! @@ -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): @@ -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. @@ -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 @@ -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 @@ -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)) diff --git a/lib/python/picongpu/extra/input/preparingInsightData_example.ipynb b/lib/python/picongpu/extra/input/preparingInsightData_example.ipynb index 53b3e6ce6e..dc1ef390e6 100644 --- a/lib/python/picongpu/extra/input/preparingInsightData_example.ipynb +++ b/lib/python/picongpu/extra/input/preparingInsightData_example.ipynb @@ -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." ] }, { @@ -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)" ] }, {