Hdf5 laser profile

CMakeLists.txt

@@ -293,6 +293,7 @@ set(PARATAXIS_VERBOSE "1" CACHE STRING "Set verbosity level for PARATAXIS")
 add_definitions(-DPARATAXIS_VERBOSE_LVL=${PARATAXIS_VERBOSE})
 option(PARATAXIS_CHECK_PHOTON_CT "Checks if the number of photons generated can fit into the memory" ON)
 option(PARATAXIS_FORCE_SINGLE_SCATTERING "Make particles scatter only once" OFF)
+option(PARATAXIS_WEIGHTED_PHOTONS "Use weighted photons" OFF)
 set(PARATAXIS_NVPROF_START_TS 200 CACHE STRING "Start profiling at this timestep")
 set(PARATAXIS_NVPROF_NUM_TS 0 CACHE STRING "Number of timesteps to profile. 0 = Disabled") 
 if(PARATAXIS_CHECK_PHOTON_CT)
@@ -305,6 +306,11 @@ if(PARATAXIS_FORCE_SINGLE_SCATTERING)
 else()
     add_definitions(-DPARATAXIS_FORCE_SINGLE_SCATTERING=0)
 endif()
+if(PARATAXIS_WEIGHTED_PHOTONS)
+    add_definitions(-DPARATAXIS_WEIGHTED_PHOTONS=1)
+else()
+    add_definitions(-DPARATAXIS_WEIGHTED_PHOTONS=0)
+endif()
 add_definitions(-DPARATAXIS_NVPROF_START_TS=${PARATAXIS_NVPROF_START_TS} -DPARATAXIS_NVPROF_NUM_TS=${PARATAXIS_NVPROF_NUM_TS})
 
 ################################################################################

examples/hdf5Test/cmakeFlags

@@ -23,9 +23,9 @@
 #   - start with zero index
 #   - increase by 1, no gaps
 
-flags[0]="-DPARATAXIS_VERBOSE=13 -DPARAM_OVERWRITES:LIST=-DPARAM_PARTICLEHANDLER=AddWaveParticles;-DPARAM_FIELD_MANIPULATOR=Random"
-flags[1]="-DPARATAXIS_VERBOSE=13 -DPARAM_OVERWRITES:LIST=-DPARAM_PARTICLEHANDLER=CountParticles;-DPARAM_FIELD_MANIPULATOR=Random"
-flags[2]="-DPARATAXIS_VERBOSE=255 -DPARAM_OVERWRITES:LIST=-DPARAM_PARTICLEHANDLER=CountParticles;-DPARAM_FIELD_MANIPULATOR=HDF5Interpolator"
+flags[0]="-DPARATAXIS_VERBOSE=13 -DPARAM_OVERWRITES:LIST=-DPARAM_PARTICLEHANDLER=AddWaveParticles;-DPARAM_FIELD_MANIPULATOR=Random;-DPARAM_LASER_PHOTON_CT=Const;-DPARAM_LASER_DISTRIBUTION=EqualToPhotons"
+flags[1]="-DPARATAXIS_VERBOSE=13 -DPARAM_OVERWRITES:LIST=-DPARAM_PARTICLEHANDLER=CountParticles;-DPARAM_FIELD_MANIPULATOR=Random;-DPARAM_LASER_PHOTON_CT=Const;-DPARAM_LASER_DISTRIBUTION=EqualToPhotons"
+flags[2]="-DPARATAXIS_VERBOSE=13 -DPARAM_OVERWRITES:LIST=-DPARAM_PARTICLEHANDLER=CountParticles;-DPARAM_FIELD_MANIPULATOR=HDF5Interpolator;-DPARAM_LASER_PHOTON_CT=HDF5Profile;-DPARAM_LASER_DISTRIBUTION=Const -DPARATAXIS_WEIGHTED_PHOTONS=ON"
 
 
 ################################################################################

examples/hdf5Test/documentation.yml

@@ -60,8 +60,15 @@ tests:
     cmakeFlag: 2
     cfgFile:   0001gpusHDF5Field.cfg
     pre-run: &createHDF5Field
-        - mkdir "${TEST_INSTALL_PATH}/testData/$TEST_NAME"
-        - python3 "${TEST_INSTALL_PATH}/testData/createHDF5Filesets.py" --folder "${TEST_INSTALL_PATH}/testData/$TEST_NAME" --size ${TEST_GRID_SIZE} --dt "2e-16" --time "1.5e-15"
+        - |
+            TESTDATA_PATH="${TEST_INSTALL_PATH}/testData"
+            TESTDATA_OUT_PATH="${TESTDATA_PATH}/$TEST_NAME"
+            mkdir "${TESTDATA_OUT_PATH}"
+            python3 "${TESTDATA_PATH}/createHDF5Filesets.py" --createDensity --folder "${TESTDATA_OUT_PATH}" --size ${TEST_GRID_SIZE} --cellSize "4e-9" "4e-9" "4e-9" --dt "2e-16" --time "1.5e-15"
+            python3 "${TESTDATA_PATH}/createHDF5Filesets.py" --createLaser --folder "${TESTDATA_OUT_PATH}" --size ${TEST_GRID_SIZE} --cellSize "4e-9" "4e-9" "4e-9" --dt "1.33e-17" --time "1.5e-15"
+            # Check for openPMD conformance
+            echo "Validating generated test files"
+            find "${TESTDATA_OUT_PATH}" -type f -name "*.h5" -print0 | xargs -0 -n 1 -r python "${TEST_INSTALL_PATH}/buildSystem/checkOpenPMD_h5.py" -i 
     post-run: &validateHDF5Field
         - *validateOpenPMD
         - python3 "${TEST_INSTALL_PATH}/testData/validateCheckpoint.py"

examples/hdf5Test/include/simulation_defines/param/laserConfig.param

@@ -0,0 +1,118 @@
+/**
+ * Copyright 2015-2016 Alexander Grund
+ *
+ * This file is part of ParaTAXIS.
+ *
+ * ParaTAXIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ParaTAXIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with ParaTAXIS.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <math/vector/Float.hpp>
+
+namespace parataxis {
+namespace laserConfig {
+
+    /* Initialization for the 'laser' */
+
+    /** Specifies the intensity of the laser (number of photons at a given position and time) */
+    namespace photonCount {
+        /* Possible values: */
+        /** Same number of particles at all points */
+        struct Const{
+            /** Value to return */
+            static constexpr float_X numPhotons = 10;
+        };
+        /** Load from HDF5 file in OpenPMD format. Requires runtime parameters */
+        struct HDF5Profile;
+
+        /** Value used in the simulation */
+        using UsedValue = PARAM_LASER_PHOTON_CT;
+
+    }  // namespace distribution
+
+    /** Specifies the distribution (number of particles at a given position and time) */
+    namespace distribution {
+        /* Possible values: */
+        /** Same number of particles at all points */
+        struct Const{
+            /** Value to return */
+            static constexpr uint32_t numParts = 10;
+        };
+
+        /** Number of Particles equals number of photons (rounded) */
+        struct EqualToPhotons;
+
+        /** Value used in the simulation */
+        using UsedValue = PARAM_LASER_DISTRIBUTION;
+
+    }  // namespace distribution
+
+    /** Specifies the in-cell position of the particles */
+    namespace position {
+
+        /* Possible values: */
+        /** Same position at all points */
+        struct Const{
+            /** In-Cell position [0,1) in x-direction */
+            static constexpr float_X X = 0.001;
+            /** In-Cell position [0,1) in y-direction */
+            static constexpr float_X Y = 0.5;
+            /** In-Cell position [0,1) in z-direction */
+            static constexpr float_X Z = 0.5;
+        };
+        /** Positions the particles evenly in the whole (possible) cell area) */
+        struct EvenDistance;
+        /** Positions the particles randomly */
+        struct Random;
+
+        /** Value used in the simulation */
+        using UsedValue = Random;
+
+    }  // namespace position
+
+    namespace phase {
+        /* Possible values: */
+        /** Wave front perpendicular to propagation, place-invariant phase */
+        struct PlaneWave;
+
+        /** Value used in the simulation */
+        using UsedValue = PlaneWave;
+    }  // namespace phase
+
+    namespace direction {
+        /* Possible values: */
+        /** Same direction at all points */
+        struct Const{
+            /** direction in x-direction */
+            static constexpr float_X X = 1;
+            /** direction in y-direction */
+            static constexpr float_X Y = 0;
+            /** direction in z-direction */
+            static constexpr float_X Z = 0;
+        };
+
+        /** Value used in the simulation */
+        using UsedValue = Const;
+    }  // namespace direction
+
+
+    /** Laser pulse length in s */
+    constexpr float_X PULSE_LENGTH = 100e-15;
+    /** Entry-Direction of the laser pulse: 0=Left->Right, 1=Top->Bottom, 2=Front->Back
+     *  Ignored for 2D, there it is always Left->Right */
+    constexpr int32_t DIRECTION = 0;
+
+}  // namespace laserConfig
+}  // namespace parataxis

examples/hdf5Test/include/simulation_defines/param/speciesDefinition.param

@@ -30,12 +30,20 @@ namespace parataxis {
 
     /*########################### define particle attributes #####################*/
 
+    value_identifier(float_X, amplitudePhotons, 1);
+
     /* describe attributes of a particle*/
     typedef PMacc::MakeSeq<
             position<position_pic>,
             direction,
+#if PARATAXIS_WEIGHTED_PHOTONS
+            amplitude<amplitudePhotons>,
+#endif
             //particleId,
-            creationTime,
+            //creationTime,
+#if PARATAXIS_FORCE_SINGLE_SCATTERING
+            wasScattered,
+#endif
             startPhase
     >::type DefaultParticleAttributes;
 
@@ -51,7 +59,6 @@ namespace parataxis {
     value_identifier(float_X, wavelengthPhotons, 0.1 * 1e-9);
     /* energy for photons in J */
     value_identifier(float_X, energyPhotons, 8 * UNITCONV_keV_to_Joule);
-    value_identifier(float_X, amplitudePhotons, 1);
     value_identifier(float_X, massPhotons, 0);
     value_identifier(float_X, chargePhotons, 0);
 
@@ -62,7 +69,9 @@ namespace parataxis {
             particleDetector<detector::PhotonDetector>,
             wavelength<wavelengthPhotons>,
             //energy<energyPhotons>,
+#if !PARATAXIS_WEIGHTED_PHOTONS
             amplitude<amplitudePhotons>,
+#endif
             mass<massPhotons>,
             charge<chargePhotons>
     > ParticleFlagsPhotons;
@@ -78,3 +87,4 @@ namespace parataxis {
     > PIC_Photons;
 
 }  // namespace parataxis
+

examples/hdf5Test/submit/0001gpusHDF5Field.cfg.in

@@ -35,8 +35,8 @@ TBG_gpu_x=1
 TBG_gpu_y=1
 TBG_gpu_z=1
 
-TBG_gridSize="-g 32 64 64"
-TBG_steps="-s 100"
+TBG_gridSize="-g 32 64 96"
+TBG_steps="-s 90"
 
 TBG_program="@PROJECT_NAME@"
 
@@ -45,14 +45,16 @@ TBG_program="@PROJECT_NAME@"
 #################################
 
 TBG_HDF5_FILE="${TEST_INSTALL_PATH}/testData/$TEST_NAME/field"
+TBG_LASER_FILE="${TEST_INSTALL_PATH}/testData/$TEST_NAME/laserProfile"
 
 TBG_plugins="--detSize 1400 900                           \
              !TBG_globalSeed                              \
              --p_position.period 0                        \
              --DensityField_printSlice.period 0           \
              --DensityField.hdf5File \"!TBG_HDF5_FILE\"   \
+             --laserSrc.p.hdf5File \"!TBG_LASER_FILE\"    \
              --PhotonDetector_print.period 0              \
-             --checkpoints 50                             \
+             --checkpoints 30                             \
              "
 
 #################################

examples/hdf5Test/submit/0008gpusHDF5Field.cfg.in

@@ -18,8 +18,8 @@ TBG_gpu_x=2
 TBG_gpu_y=2
 TBG_gpu_z=2
 
-TBG_gridSize="-g 64 128 128"
-TBG_steps="-s 100"
+TBG_gridSize="-g 32 64 96"
+TBG_steps="-s 90"
 
 TBG_program="@PROJECT_NAME@"
 
@@ -28,14 +28,16 @@ TBG_program="@PROJECT_NAME@"
 #################################
 
 TBG_HDF5_FILE="${TEST_INSTALL_PATH}/testData/$TEST_NAME/field"
+TBG_LASER_FILE="${TEST_INSTALL_PATH}/testData/$TEST_NAME/laserProfile"
 
 TBG_plugins="--detSize 1400 900                           \
              !TBG_globalSeed                              \
              --p_position.period 0                        \
              --DensityField_printSlice.period 0           \
              --DensityField.hdf5File \"!TBG_HDF5_FILE\"   \
+             --laserSrc.p.hdf5File \"!TBG_LASER_FILE\"    \
              --PhotonDetector_print.period 0              \
-             --checkpoints 50                             \
+             --checkpoints 30                             \
              "
 
 #################################

examples/hdf5Test/testData/createHDF5Filesets.py

@@ -28,7 +28,7 @@ def writeOpenPMDHeader(f, basePath):
     attr = f.attrs
     # Required attributes
     attr["openPMD"] = np.string_("1.0.0")
-    attr["openPMDextension"] = 0
+    attr["openPMDextension"] = np.uint32(0)
     attr["basePath"] = np.string_("/data/%T/")
     attr["meshesPath"] = np.string_("meshes/")
     attr["particlesPath"] = np.string_("particles/")
@@ -50,7 +50,7 @@ def createBaseGroup(f, timestep, dt):
     baseGroup.attrs["timeUnitSI"] = 1.0
     return baseGroup
 
-def createDensityFileset(basePath, shape, numTimesteps, dt, func):
+def createDensityFileset(basePath, shape, cellSizes, numTimesteps, dt, func):
     """Create a set of HDF5 files containing an electron_density field of the given 3D shape for the 
        period of numTimesteps of length dt with values of func(x, y, z, t)"""
     for timestep in range(numTimesteps):
@@ -59,41 +59,94 @@ def createDensityFileset(basePath, shape, numTimesteps, dt, func):
             writeOpenPMDHeader(f, basePath)
             baseGroup = createBaseGroup(f, timestep, dt)
             meshes = baseGroup.create_group(f.attrs["meshesPath"])
+
             ds = meshes.create_dataset("electron_density", data = data)
-            ds.attrs["dataOrder"] = np.string_("C")
             ds.attrs["axisLabels"] = np.array([np.string_("z"), np.string_("y"), np.string_("x")])
+            ds.attrs["dataOrder"] = np.string_("C")
+            ds.attrs["geometry"] = np.string_("cartesian")
+            ds.attrs["unitDimension"] = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+            ds.attrs["gridGlobalOffset"] = np.array([0.0, 0.0, 0.0])
+            gridUnitSI = 1.
+            ds.attrs["gridSpacing"] = cellSizes / gridUnitSI
+            ds.attrs["gridUnitSI"] = gridUnitSI
+            ds.attrs["timeOffset"] = 0.0
+            ds.attrs["position"] = np.array([0.0, 0.0, 0.0])
+            ds.attrs["unitSI"] = 1.0
+
 
-def createPhotonFileset(basePath, shape, numTimesteps, dt, func):
+def createPhotonFileset(basePath, shape, cellSizes, numTimesteps, dt, func):
     """Create a set of HDF5 files containing the photon count in the given 2D shape for the 
        period of numTimesteps of length dt with values of func(x, y, t)"""
     for timestep in range(numTimesteps):
         data = np.fromfunction(np.vectorize(lambda y,x: func(x, y, timestep * dt)), shape)
         with h5py.File(basePath + "_" + str(timestep) + ".h5", "w") as f:
             writeOpenPMDHeader(f, basePath)
             baseGroup = createBaseGroup(f, timestep, dt)
+
             numPhotons = baseGroup.create_group(f.attrs["meshesPath"]).create_group("Nph")
+            numPhotons.attrs["axisLabels"] = np.array([np.string_("y"), np.string_("x")])
+            numPhotons.attrs["dataOrder"] = np.string_("C")
+            numPhotons.attrs["geometry"] = np.string_("cartesian")
+            numPhotons.attrs["gridGlobalOffset"] = np.array([0.0, 0.0])
+            # Choose an arbitrary unit as additional test
+            gridUnitSI = 2.
+            numPhotons.attrs["gridSpacing"] = cellSizes / gridUnitSI
+            numPhotons.attrs["gridUnitSI"] = gridUnitSI
+            numPhotons.attrs["timeOffset"] = 0.0
+            numPhotons.attrs["unitDimension"] = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+
+            # x-polarized
             ds = numPhotons.create_dataset("x", data = data)
-            ds.attrs["dataOrder"] = np.string_("C")
-            ds.attrs["axisLabels"] = np.array([np.string_("y"), np.string_("x")])
+            ds.attrs["position"] = np.array([0.0, 0.0])
+            ds.attrs["unitSI"] = 1.0
+
+            # y-polarized (nothing yet)
             ds = numPhotons.create_dataset("y", data = np.zeros(shape))
-            ds.attrs["dataOrder"] = np.string_("C")
-            ds.attrs["axisLabels"] = np.array([np.string_("y"), np.string_("x")])
+            ds.attrs["position"] = np.array([0.0, 0.0])
+            ds.attrs["unitSI"] = 1.0
 
 # Parse the command line.
 parser = argparse.ArgumentParser(description="Create an HDF5 file set for the simulation", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument('--createDensity', action='store_true', help='Create electron density')
+parser.add_argument('--createLaser', action='store_true', help='Create XRay photon distribution')
 parser.add_argument('-f', '--folder', help='Output path')
 parser.add_argument('-s', '--size', type=int, nargs=3, help='Size of the simulation (x y z)')
+parser.add_argument('-c', '--cellSize', type=float, nargs=3, help='Size of the simulation cells (x y z)')
 parser.add_argument('--dt', type=float, help='Timestep used for HDF5 records')
 parser.add_argument('-t', '--time', type=float, help='Total time for HDF5 records')
 options = parser.parse_args()
 
 numTimesteps = int(math.ceil(options.time / options.dt))
 size = np.flipud(np.array(options.size)) # z y x (fastest varying index last)
-size2D = np.array([size[1], size[0]])
-# Create time varying circle
-createDensityFileset(os.path.join(options.folder, "field"), size, numTimesteps, options.dt,
-                     lambda x,y,z,t: t * 0.65/options.dt if np.linalg.norm((y,z) - size2D/2) < 10 else 0.)
-# Create time varying circle
-createPhotonFileset(os.path.join(options.folder, "laserProfile"), size2D, numTimesteps, options.dt,
-                     lambda x,y,t: t/options.dt * max(15 - np.linalg.norm((x,y) - size2D/2), 0))
+size2D = np.array([options.size[1], options.size[2]]) # Y=Y Z=X
+cellSize = np.flipud(np.array(options.cellSize)) # z y x (fastest varying index last)
+cellSize2D = np.array([options.cellSize[1], options.cellSize[2]]) # Y=Y Z=X
+
+def calcDensityField(x, y, z, t, dt):
+    # Don't scatter in first few cells to allow laser test
+    if x < 5:
+        return 0.
+    if np.linalg.norm((y,z) - size2D/2) < 10:
+        return t * 0.65/options.dt
+    else:
+        return 0.
+
+def getPhotonCount(idxX, idxY, timestep):
+    """Return the number of photons expected in the given cell and timestep"""
+    result = max(idxX - 4, 0) + max(idxY - 2, 0) * 3 + max(timestep - 4, 0)
+    return float(result)
+
+if options.createDensity:
+    # Create time varying circle
+    createDensityFileset(os.path.join(options.folder, "field"), size, cellSize, numTimesteps, options.dt,
+                         lambda x,y,z,t: calcDensityField(x, y, z, t, options.dt) )
+
+if options.createLaser:
+    # Here we can make the hdf5 cells bigger than the sim cells
+    xRatio = 1
+    yRatio = 1
+    tRatio = 1
+    # the function is defined in terms of simulation units, so multiply our cells and timestep and then scale with area and time ratio
+    photonFunc = lambda x,y,t: getPhotonCount(x*xRatio, y*yRatio, int(t/options.dt)*tRatio) * xRatio*yRatio*tRatio
+    createPhotonFileset(os.path.join(options.folder, "laserProfile"), size2D, cellSize2D * np.array([yRatio, xRatio]), numTimesteps, options.dt * tRatio, photonFunc)