High performance Amplified Spontaneous Emission on GPU
HASEonGPU is a scientific project. If you present and/or publish scientific results that used HASEonGPU, you should set this as a reference.
HASEonGPU is licensed under the GPLv3+. Please refer to our LICENSE.md
- 
Software:
- make
 - cmake 2.8.10
 - gcc 4.4.1
 - cuda 5.5
 
 - 
Optional:
- octave / matlab
 - paraview
 - OpenMPI 1.8 or compatible
 
 - 
Hardware:
- Nvidia device >= Compute capability 2.0 (at least fermi generation)
 
 
- clone the repository: 
git clone https://github.com/computationalradiationphysics/haseongpu.git - create the build directory: 
mkdir haseongpu/build - go to build directory: 
cd haseongpu/build - create Makefile 
cmake .. - build project : 
make 
- MATLAB compatible interface
 - C-Application interface
 
A small example for Phi ASE calculation with a pumped crystal. The simulation can be started by the following:
- follow the compile instructions above
 - change path "cd example/matlab_example/"
 - run : 
matlab laserPumpCladdingExample - watch progress
 - take a look at the results (*.vtk) with paraview
 
- follow the compile instructions above
 - change path "cd example/c_example/"
 - run : 
./bin/calcPhiASE --input-path=./input/cylindrical --min-rays=10000 - watch progress
 - take a look at the results in the output directory
 
- Add calcPhiASE.m to your matlab path
 - Call calcPhiASE from your matlab script like following:
 
[phiASE, MSE, nRays] = calcPhiASE(
                                    points,  
                                    trianglePointIndices,           
                                    betaCells,  
                                    betaVolume,  
                                    claddingCellTypes,  
                                    claddingNumber,  
                                    claddingAbsorption,  
                                    useReflections,  
                                    refractiveIndices,  
                                    reflectivities,  
                                    triangleNormalsX,  
                                    triangleNormalsY,  
                                    triangleNeighbors,  
                                    triangleSurfaces,  
                                    triangleCenterX,  
                                    triangleCenterY,  
                                    triangleNormalPoint,  
                                    forbiddenEdge,  
                                    minRaysPerSample,  
                                    maxRaysPerSample,  
                                    mseThreshold,  
                                    repetitions,  
                                    nTot,  
                                    thickness,  
                                    laserParameter,  
                                    crystal,  
                                    numberOfLevels,  
                                    deviceMode,
                                    parallelMode,  
                                    maxGPUs,  
                                    nPerNode
                                );  - 
The returned values are represented as two-dimensional matrices in which columns are slice indices(levels) and rows are point indices. The value for the ith point and jth slice can then be optained by MATLAB with:
value = values(i,j);
 
In the following all arguments of the MATLAB call are described. You will find on each point a head with datatype (an array when in brackets []), the size of the array and to which set of numbers the array belongs.
- 
points [float], in {0, ..., numberOfPoints}, size = numberOfPoints
The coordinates of the triangle vertices. All x coordinates followed by all of the y coordinates of the triangle vertices
structure: [x_1, x_2, ... x_n, y_1, y_2, ... y_n] (n == numberOfPoints) - 
trianglePointIndices [int] in {0, ..., numberOfPoints}, size = numberOfTriangles * 3
Contains the indices to access the "points" datastructure (each triangle has 3 points as vertices). Each entry is an index from 0 to numberOfPoints, corresponding to the positions of a vertex in "points". Structure as follows:
[ triangle1A, triangle2A, ... triangleNA, triangle1B, triangle2B, ... triangleNB, triangle1C, ... ]
i.e. for triangles with vertices A,B,C there are all the indices of the A-vertices, followed by all the B and C vertices. - 
betaCells [float]
Stimulus in the sample points. - 
betaVolume [float], size = numberOfTriangles * numberOfLevels - 1
Stimulus in the volume (prisms). Beta values for all prisms ordered accordingly to the prismIDs: prismID = triangleID + layer * numberOfTriangles. Therefore, all betaValues for a layer are grouped together - 
claddingCellTypes [int], size = numberOfTriangles
Sets cladding index for triangles {0,1,2,...} - 
claddingNumber unsigned, size = 1
Set which cladding to use - 
claddingAbsorption float, size = 1
Absorption coefficient of cladding - 
useReflections bool, size = 1
Switch to activate reflections. - 
refractiveIndices [float], size = 4
Describes the refractive indices of the active gain medium top and bottom planes. It is structured as follows: {bottomInside, bottomOutside, topInside, topOutside} bottomInside = topInside (because it is the same medium) - 
reflectivities [float], in {0, ...,1}, size = 2 * numberOfTriangles
Defines the reflectivities of prism planes. First the reflectivities of bottom plane and then the reflectivities of top plane. Both it ordered by TriangleID. - 
triangleNormalsX [float], size = numberOfTriangles * 3
The x coordinate of the normal vectors for each triangle edge. It is ordered as follows:
[ triangle1_edge0, triangle2_edge0, ... triangleN_edge0, triangle1_edge1, triangle2_edge1, ... ]
i.e. all first edges of each triangle, followed by all second edges of each triangle and so on. - 
triangleNormalsY [float], size = numberOfTriangles * 3
The y coordinate of the normal vectors for each triangle edge. It is ordered as follows:
[ triangle1_edge0, triangle2_edge0, ... triangleN_edge0, triangle1_edge1, triangle2_edge1, ... ]
i.e. all first edges of each triangle, followed by all second edges of each triangle and so on. - 
triangleNeighbors [int], in {-1,0,1,2,4}, size = 5
Describes the neighnor relation of triangles to each other. Each entry corresponds to a triangleID (see "triangles") which is adjacent to the current triangle and edge. Structure is similar to "forbidden":
[ triangle1_edge0, triangle2_edge0, ... triangleN_edge0, triangle1_edge1, triangle2_edge1, ... ] - 
triangleSurfaces [float], size = numberOfTriangles
The sizes of the surfaces of each triangle, ordered by the triangleID. - 
triangleCenterX [float], size = numberOfTriangles
The x coordinates of the center points for each triangle ordered by TriangleID. - 
triangleCenterY [float], size = numberOfTriangles
The y coordinates of the center points for each triangle ordered by TriangleID. - 
triangleNormalPoint [unsigned], in {0, ..., numberOfPoints}, size = numberOfTriangles * 3
Contains indices to the point where the triangleNormalVectors start. For each Triangle 3 points (3 edges) are stored in this list. Indices point to locations in "points" (i.e. normal vectors start at triangle vertices!)l Structure is VERY similar to triangles: [ triangle1_p0, triangle2_p0, ... triangleN_p0, triangle1_p1, triangle2_p1, ... ] - 
forbiddenEdge [int], in {-1,0,1,2,4}, size = 5
Describes the relation of edge indices of adjacent triangles -1 means, there is no adjacent triangle to that edge 0,1,2 describes the index of the edge as seen from the ADJACENT triangle Order of data is similar to normalVec: [ triangle1_edge0, triangle2_edge0, ... triangleN_edge0, triangle1_edge1, triangle2_edge1, ... ] i.e. all first edges of each triangle, followed by all second edges of each triangle and so on. - 
minRaysPerSample unsigned, size = 1
Minimal number of rays for adaptive sampling - 
maxRaysPerSample unsigned, size = 1
Maximal number of rays for adaptive sampling - 
mseThreshold float, size = 1
Sets the maximal MSE of the ASE value. If a sample-point does not reach this MSE-threshold, the number of rays per sample-point will be increased upto maxRaysPerSample or resampled with repetitive sampling. - 
repetitions unsigned, size = 1
Sets the number of maximal repetitions when the mseThreshold was not reached. - 
nTot float, size = 1
Doping of the active gain medium - 
thickness float, size = 1
Thickness of one prism level of the mesh. - 
laserParameter [float]
Is a structure for the laser parameters (intensities sigma, wavelength lambda) s_ems corresponds to l_ems and s_abs to l_abs struct(s_abs, VALUES, s_ems, VALUES, l_abs, VALUES, l_ems, VALUES) - 
crystal [float]
Is a structure for the crystal parameters crystal.tfluo describes the crystalFluorescence of the active gain medium. - 
numberOfLevels unsigned, size = 1
Total number of levels of the mesh. Thus the total thickness of the mesh is thickness * numberOfLevels! - 
parallelMode
+ 'mpi' use mpi to distribute workload (use nPerNode) + 'threaded' use pthreads to distribute workload locally - 
deviceMode + 'cpu' use cpu algorithm (does not have all features) + 'gpu' use gpu algorithm
 - 
maxGpus unsigned, size = 1
Maximal number of GPUs for threaded case - 
nPerNode
Number of devices per mpi-node 
- 
Command:
./bin/calcPhiASE [OPTIONS]
 - 
Options:
 
--input-path 
  Path to the experiment location.
  This folder contains several .txt files usually
  generated by an matlab script. The content of this
  .txt files contains all experiment data you need
  to run one experiment.
    
--output-path
  Path to a writable location. Is used to write
  input and output for matlab script.
--parallel-mode=[|threaded|mpi]  
  Defines the method of parallelization to start the
  simulation with. Mode "threaded" uses pthreads on a single
  node. Mode "mpi" is a parallel mpi
  implementation for clusters. Note, that this parameter
  is currently only available when using `--device-mode=gpu`
--device-mode=[cpu|gpu]  
  Defines on which hardware the simulation will run.
  Mode "cpu" is the original
  algorithm based on single core cpu.
  Mode "gpu" uses nVIDIA CUDA GPUs, that can be parallelized either
  with Pthreads or MPI.
--min-rays=  
  Sets the minimum number of rays per sample-point in the
  crystal structure.
--max-rays=  
  Sets the maximal number of rays per sample-point. The number
  of rays per sample-point will vary between minimum and
  maximum number of rays in dependance of a MSE-Threshold.
  (see --mse-threshold)
--ngpus=  
  Set the number of gpus to use. "mpi" parallel-mode should set this
  to 1 and "threaded" to the maximum number
  of GPUs on the node. If you don't set it, it will
  be set to the maximum automatically.
--min-sample-i=  
  Index of the first sample point (normally 0).
--max-sample-i=  
  Index of the last sample point (numberOfSample - 1).
      
--verbosity=  
  Add the following for different verbosity levels:
  0  : quiet
  1  : error
  2  : warning
  4  : info
  8  : statistics
  16 : debug
  32 : progress-bar
  Levels the verbosity level is interpreted as a bitmask and 
  can be composed of different levels.
--reflection  
  Use reflection on upper and lower plane of gain
  medium. Maximal number of reflections will be
  calculated
--mse-threshold=  
  Algorithm tries to stay under this threshold
  by adaptive and repetitive sampling.
--repetitions=  
  Number of repetitions, that will be done
  when mse-threshold was not met.
--spectral-resolution= 
  Resolution of absorption and emission spectrum to which the
  input spectrum will be interpolated linear.  Interpolation is
  used to distribute spectrum values equidistant over the
  wavelength.  Omitting this option or setting a to small
  resolutionwill set the lambda resolution to the maximum number
  of absorption or emission values.
- 
4 GPUs, 10K to 100K Rays, 4 Repetitions
./bin/calcPhiASE --input-path=/input/ --output-path=/tmp/ --parallel-mode=threaded --min-rays=10000 --max-rays=100000 --reflection --repetitions=4 --ngpus=4 --min-sample-i=0 --max-sample-i=1234 --mse-threshold=0.05 - 
MPI with 4 GPUs per node
mpiexec -npernode 4 ./bin/calcPhiASE --input-path=/input/ --output-path=/tmp/ --parallel-mode=mpi --min-rays=10000 --max-rays=100000 --reflection --repetitions=4 --ngpus=1 --min-sample-i=0 --max-sample-i=1234 --mse-threshold=0.05 
- Dr. Michael Bussmann
 - Dr. Daniel Albach
 
- Erik Zenker
 - Carlchristian Eckert
 
- Marius Melzer
 - Frank Liebold
 - Daniel Höhne
 
CMakeLists.txtcmake file to generate a Makefilesrc/folder containing all the source code that is not a headersrc/map_rays_to_prisms.cuCUDA code to generate a schedule of which ray will be launched from which prism of the gain mediumsrc/calcPhiASE.mMATLAB adapter scriptsrc/logging.cucreates nicely readable output based on log-levelssrc/importance_sampling.cuCUDA parallelized importance samplingsrc/geometry.cubasic 3D geometry calculationssrc/interpolation.cuinterpolation functions for wavelengths of polychromatic laser pulsessrc/reflection.cuCUDA functions to calculate reflections inside the gain mediumsrc/progressbar.cuprogressbar for the command linesrc/parser.cuparsing of command line arguments and textfilessrc/for_loops_clad.cuold CPU code for ASE calculationsrc/calc_phi_ase_mpi.ccMPI workload distribution. Code for Master and Slavessrc/write_to_file.cuwriting formatted data to a filesrc/ray_histogram.cuprint a histogram of the adaptive ray count to command linesrc/calc_phi_ase_threaded.cupthreads workload distributionsrc/mt19937ar.cuCPU code for Mersenne Twister PRNG used by for_loops_clad.cusrc/write_to_vtk.cugenerate VTK-files from simulation results (deprecated)src/propagate_ray.cuCUDA code to propagate a single ray through the prism mesh structuresrc/mesh.cuclass that holds the information and all parameters about the gain medium meshsrc/cuda_utils.cuutility functions (getFreeDevice)src/calc_sample_gain_sum.cuCUDA code to calculate all the rays for a single sample pointsrc/calc_phi_ase.cuCUDA code to calculate ASE for all the sample pointssrc/main.cumain entry filesrc/write_matlab_output.cugenerate MATLAB-readable matrixes from the simulation dataexample/folder that contains executable examplesexample/c_example/folder that contains the commandline-exampleexample/c_example/input/folder that contains input for 2 different experiments -example/c_example/input/cylindrical/folder that contains data for the cylindrical gain medium. For details on the files, see detailed information above (Input argument description) -example/c_example/input/cuboid/example input with a cuboid gain medium. contents similar to cylindrical example.example/c_example/output/folder to gather the outputexample/matlab_example/folder that contains the input data for the matlab exampleexample/matlab_example/lambda_e.txtemission wavelengthsexample/matlab_example/sigma_e.txtemission crosssectionexample/matlab_example/pt.matsampling points and delaynay-triangles of the gain mediumexample/matlab_example/set_variables.mgenerate information about the meshexample/matlab_example/vtk_wedge.mgenerate a VTK file from the meshexample/matlab_example/laserPumpCladdingExample.mexperimental setup. Run this file to see the progress of a whole experimentexample/matlab_example/sve.matexample/matlab_example/sigma_a.txtabsorption crosssectionexample/matlab_example/gain.mcalculate gain distribution inside the gain mediumexample/matlab_example/beta_int3.mutility function to calculate gain distributionexample/matlab_example/extract_gain_map.mcalculate the gain for the sample point used in the actual measurementexample/matlab_example/beta_int.mutility function to calculate gain distributionexample/matlab_example/lambda_a.txtabsorption wavelengthsinclude/folder containing all the header source codeinclude/calc_phi_ase_mpi.hppheader for calc_phi_ase_mpi.cuinclude/mesh.hppheader for mesh.cuinclude/importance_sampling.hppheader for importance_sampling.cuinclude/ray_histogram.hppheader for ray_histogram.cuinclude/for_loops_clad.hppheader for for_loops_clad.cuinclude/mt19937ar.hppheader for mt19937ar.cuinclude/calc_phi_ase.hppheader for calc_phi_ase.cuinclude/write_matlab_output.hppheader for write_matlab_output.cuinclude/cuda_utils.hppheader for cuda_utils.cuinclude/logging.hppheader for logging.cuinclude/cudachecks.hppMacros to check the success state of CUDA callsinclude/reflection.hppheader for reflection.cuinclude/parser.hppheader for parser.cuinclude/map_rays_to_prisms.hppheader for map_rays_to_prisms.cuinclude/calc_phi_ase_threaded.hppheader for calc_phi_ase_threaded.cuinclude/thrust_device_vector_nowarn.hppwrapper to switch off compiler warning that is produced by 3rd party library (CUDA Thrust)include/propagate_ray.hppheader for propagate_ray.cuinclude/thrust_host_vector_nowarn.hppwrapper to switch off compiler warning that is produced by 3rd party library (CUDA Thrust)include/calc_sample_gain_sum.hppheader for calc_sample_gain_sum.cuinclude/interpolation.hppheader for interpolation.cuinclude/version.hppversion information for HASEonGPUinclude/geometry.hppheader for geometry.cuinclude/write_to_file.hppheader for write_to_file.cuinclude/types.hpptype definitions for HASEonGPUinclude/progressbar.hppheader for progressbar.cuinclude/nan_fix.hppwrapper to allow usage ofisnan()in a templateinclude/write_to_vtk.hppheader for write_to_vtk.cuLICENSE.mdadditional licensing informationREADME.mdthis README fileREFERENCE.mdReferencing informationCOPYINGFull License informationutils/folder that contains utility filesutils/cmake/utils/cmake/modules/3rd Party CMAKE module that was modified to circumvent a bug where the NVCC linker would crash on unknown arguments