codebase / stimlib

Commit 4d67ff4e34005520e69b9cb28a0f603f5abd4bde

Authored by David Mayerich
1 parent 3b012a80

started scalar field and efield classes

Showing 3 changed files with 220 additions and 0 deletions   Show diff stats
Inline Side-by-side
source/colormap.cpp renamed to cpp/colormap.cpp
  Show/Hide comments   View file @4d67ff4
optics/efield.cuh 0 → 100644
  Show/Hide comments   View file @4d67ff4
  1 +#include "rts/math/complex.h"
  2 +#include "rts/visualization/colormap.h"
  3 +
  4 +namespace rts{
  5 +
  6 +/* This class implements a discrete representation of an electromagnetic field
  7 + in 2D. The majority of this representation is done on the GPU.
  8 +*/
  9 +template<typename P>
  10 +class efield
  11 +{
  12 +private:
  13 +
  14 + bool scalar;
  15 +
  16 + //gpu pointer to the field data
  17 + rts::complex<P>* X;
  18 + rts::complex<P>* Y;
  19 + rts::complex<P>* Z;
  20 +
  21 + //resolution of the discrete field
  22 + unsigned int R[2];
  23 +
  24 +
  25 +public:
  26 +
  27 + efield(unsigned int res0, unsigned int res1, bool _scalar = false)
  28 + {
  29 + //initialize all pointers to NULL
  30 + X = Y = Z = NULL;
  31 + R[0] = res0;
  32 + R[1] = res1;
  33 +
  34 + //allocate memory on the gpu
  35 + cudaMalloc(&X, sizeof(rts::complex<P>) * R[0] * R[1]);
  36 +
  37 + if(!scalar)
  38 + {
  39 + cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
  40 + cudaMalloc(&Z, sizeof(rts::complex<P>) * R[0] * R[1]);
  41 + }
  42 + }
  43 +
  44 + //destructor
  45 + ~efield()
  46 + {
  47 + if(X != NULL) cudaFree(X);
  48 + if(Y != NULL) cudaFree(Y);
  49 + if(Z != NULL) cudaFree(Z);
  50 + }
  51 +
  52 + void render(std::string)
  53 + {
  54 +
  55 + }
  56 +
  57 +
  58 +};
  59 +
  60 +
  61 +
  62 +} //end namespace rts
... ...
visualization/scalarfield.cuh 0 → 100644
  Show/Hide comments   View file @4d67ff4
  1 +#ifndef RTS_SCALAR_SLICE
  2 +#define RTS_SCALAR_SLICE
  3 +
  4 +#include "rts/visualization/colormap.h"
  5 +#include "rts/envi/envi.h"
  6 +#include "cublas_v2.h"
  7 +#include <cuda_runtime.h>
  8 +
  9 +template<typename P>
  10 +struct scalarfield
  11 +{
  12 + //gpu pointer to the scalar slice
  13 + P* S;
  14 +
  15 + //resolution of the slice
  16 + int R[2];
  17 +
  18 + scalarfield()
  19 + {
  20 + R[0] = R[1] = 0;
  21 + S = NULL;
  22 +
  23 + //std::cout<<"Scalerslice created (default)."<<std::endl;
  24 + }
  25 +
  26 + scalarfield(int x, int y)
  27 + {
  28 + //set the resolution
  29 + R[0] = x;
  30 + R[1] = y;
  31 +
  32 + //allocate memory on the GPU
  33 + HANDLE_ERROR(cudaMalloc( (void**)&S, sizeof(P) * x * y ));
  34 + //std::cout<<"Scalerslice created."<<std::endl;
  35 + }
  36 +
  37 + ~scalarfield()
  38 + {
  39 + if(S != NULL)
  40 + HANDLE_ERROR(cudaFree(S));
  41 + S = NULL;
  42 + R[0] = R[1] = 0;
  43 +
  44 + //std::cout<<"Scalerslice destroyed."<<std::endl;
  45 + }
  46 +
  47 + void clear()
  48 + {
  49 + //this function sets the slice to zero
  50 + if(S != NULL)
  51 + //HANDLE_ERROR(cudaMalloc( (void**)&S, sizeof(P) * R[0] * R[1] ));
  52 + HANDLE_ERROR(cudaMemset(S, 0, sizeof(P) * R[0] * R[1]));
  53 + }
  54 +
  55 + void toImage(std::string filename, P vmin, P vmax, rts::colormapType cmap = rts::cmBrewer)
  56 + {
  57 + rts::gpu2image<ptype>(S, filename, R[0], R[1], vmin, vmax, cmap);
  58 + }
  59 +
  60 + void toImage(std::string filename, bool positive = true, rts::colormapType cmap = rts::cmBrewer)
  61 + {
  62 + cublasStatus_t stat;
  63 + cublasHandle_t handle;
  64 +
  65 + //create a CUBLAS handle
  66 + stat = cublasCreate(&handle);
  67 + if(stat != CUBLAS_STATUS_SUCCESS)
  68 + {
  69 + std::cout<<"CUBLAS Error: initialization failed"<<std::endl;
  70 + exit(1);
  71 + }
  72 +
  73 + //find the index of the value with maximum magnitude
  74 + int N = R[0] * R[1];
  75 + int result;
  76 + #ifdef PRECISION_SINGLE
  77 + stat = cublasIsamax(handle, N, S, 1, &result);
  78 + #elif defined PRECISION_DOUBLE
  79 + stat = cublasIdamax(handle, N, S, 1, &result);
  80 + #endif
  81 +
  82 + //adjust for 1-based indexing
  83 + result -= 1;
  84 +
  85 + if(stat != CUBLAS_STATUS_SUCCESS)
  86 + {
  87 + std::cout<<"CUBLAS Error: failure finding maximum value."<<std::endl;
  88 + exit(1);
  89 + }
  90 +
  91 +
  92 +
  93 + //retrieve the maximum value
  94 + P maxVal;
  95 + HANDLE_ERROR(cudaMemcpy(&maxVal, S + result, sizeof(ptype), cudaMemcpyDeviceToHost));
  96 +
  97 + //destroy the CUBLAS handle
  98 + cublasDestroy(handle);
  99 +
  100 + //output the image
  101 + if(positive)
  102 + toImage(filename, 0, maxVal, cmap);
  103 + else
  104 + toImage(filename, -abs(maxVal), abs(maxVal), cmap);
  105 + }
  106 +
  107 +
  108 + void toEnvi(std::string filename, ptype wavelength = 0, bool append = false)
  109 + {
  110 + std::string mode;
  111 + if(append) mode = "a";
  112 + else mode = "w";
  113 +
  114 + //open the ENVI file
  115 + EnviFile outfile(filename, mode);
  116 +
  117 + //get the scalar slice from the GPU to the CPU
  118 + int memsize = sizeof(ptype) * R[0] * R[1];
  119 + ptype* cpuData = (ptype*) malloc( memsize );
  120 + HANDLE_ERROR(cudaMemcpy( cpuData, S, memsize, cudaMemcpyDeviceToHost));
  121 +
  122 + //add a band to the ENVI file
  123 + outfile.addBand(cpuData, R[0], R[1], wavelength);
  124 +
  125 + outfile.close();
  126 +
  127 +
  128 + }
  129 +
  130 + //assignment operator
  131 + scalarfield & operator= (const scalarfield & rhs)
  132 + {
  133 + //de-allocate any existing GPU memory
  134 + if(S != NULL)
  135 + HANDLE_ERROR(cudaFree(S));
  136 +
  137 + //copy the slice resolution
  138 + R[0] = rhs.R[0];
  139 + R[1] = rhs.R[1];
  140 +
  141 + //allocate the necessary memory
  142 + HANDLE_ERROR(cudaMalloc(&S, sizeof(ptype) * R[0] * R[1]));
  143 +
  144 + //copy the slice
  145 + HANDLE_ERROR(cudaMemcpy(S, rhs.S, sizeof(ptype) * R[0] * R[1], cudaMemcpyDeviceToDevice));
  146 +
  147 +
  148 + std::cout<<"Assignment operator."<<std::endl;
  149 +
  150 + return *this;
  151 +
  152 + }
  153 +
  154 +};
  155 +
  156 +
  157 +
  158 +#endif
... ...