Skip to content

Commit

Permalink
initial port of mmc momentum transfer calculations for DCS
Browse files Browse the repository at this point in the history
  • Loading branch information
@fangq
fangq committed Apr 7, 2018
1 parent acef0cf commit 63bce2b
Show file tree
Hide file tree
Showing 4 changed files with 26 additions and 11 deletions.
22 changes: 14 additions & 8 deletions src/mcx_core.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,13 +205,13 @@ __device__ inline void savedetphoton(float n_det[],uint *detectedphoton,float ns
uint i;
for(i=0;i<gcfg->issaveseed*RAND_BUF_LEN;i++)
seeddata[baseaddr*RAND_BUF_LEN+i]=t[i]; ///< save photon seed for replay
baseaddr*=gcfg->maxmedia+2+gcfg->issaveexit*6;
baseaddr*=gcfg->maxmedia+2+gcfg->issaveexit*6+(gcfg.ismomentum*gcfg->maxmedia);
n_det[baseaddr++]=detid;
n_det[baseaddr++]=nscat;
for(i=0;i<gcfg->maxmedia;i++)
for(i=0;i<gcfg->maxmedia*(1+gcfg->ismomentum);i++)
n_det[baseaddr+i]=ppath[i]; ///< save partial pathlength to the memory
if(gcfg->issaveexit){
baseaddr+=gcfg->maxmedia;
baseaddr+=gcfg->maxmedia*(1+gcfg->ismomentum);
*((float3*)(n_det+baseaddr))=float3(p0->x,p0->y,p0->z);
baseaddr+=3;
*((float3*)(n_det+baseaddr))=float3(v->x,v->y,v->z);
Expand Down Expand Up @@ -1016,8 +1016,8 @@ kernel void mcx_main_loop(uint media[],float field[],float genergy[],uint n_seed
#endif

#ifdef SAVE_DETECTORS
ppath+=threadIdx.x*gcfg->maxmedia; // block#2: maxmedia*thread number to store the partial
if(gcfg->savedet) clearpath(ppath,gcfg->maxmedia);
ppath+=threadIdx.x*gcfg->maxmedia*(1+gcfg->ismomentum); // block#2: maxmedia*thread number to store the partial
if(gcfg->savedet) clearpath(ppath,gcfg->maxmedia*(1+gcfg->ismomentum));
#endif

gpu_rng_init(t,n_seed,idx);
Expand Down Expand Up @@ -1091,7 +1091,11 @@ kernel void mcx_main_loop(uint media[],float field[],float genergy[],uint n_seed
sincosf(theta,&stheta,&ctheta);
}
GPUDEBUG(("scat theta=%f\n",theta));

#ifdef SAVE_DETECTORS
/** accummulate momentum transfer */
if(gcfg->ismomentum)
ppath[gcfg->maxmedia+(mediaid & MED_MASK)-1]+=1.f-ctheta;
#endif
/** Update direction vector with the two random angles */
if(gcfg->is2d)
rotatevector2d(v,stheta,ctheta);
Expand Down Expand Up @@ -1541,15 +1545,15 @@ void mcx_run_simulation(Config *cfg,GPUInfo *gpu){
uint *gPseed,*gdetected;
float *gPdet,*gsrcpattern,*gfield,*genergy,*greplayw=NULL,*greplaytof=NULL,*gdebugdata=NULL;
RandType *gseeddata=NULL;
int detreclen=cfg->medianum+1+(cfg->issaveexit>0)*6;
int detreclen=cfg->medianum+1+(cfg->issaveexit>0)*6+(cfg->ismomentum*(cfg->medianum-1));
unsigned int is2d=(cfg->dim.x==1 ? 1 : (cfg->dim.y==1 ? 2 : (cfg->dim.z==1 ? 3 : 0)));
MCXParam param={cfg->steps,minstep,0,0,cfg->tend,R_C0*cfg->unitinmm,
(uint)cfg->issave2pt,(uint)cfg->isreflect,(uint)cfg->isrefint,(uint)cfg->issavedet,1.f/cfg->tstep,
p0,c0,maxidx,uint3(0,0,0),cp0,cp1,uint2(0,0),cfg->minenergy,
cfg->sradius*cfg->sradius,minstep*R_C0*cfg->unitinmm,cfg->srctype,
cfg->srcparam1,cfg->srcparam2,cfg->voidtime,cfg->maxdetphoton,
cfg->medianum-1,cfg->detnum,cfg->maxgate,0,0,cfg->reseedlimit,ABS(cfg->sradius+2.f)<EPS /*isatomic*/,
(uint)cfg->maxvoidstep,cfg->issaveseed>0,cfg->issaveexit>0,cfg->issaveref>0,
(uint)cfg->maxvoidstep,cfg->issaveseed>0,cfg->issaveexit>0,cfg->issaveref>0,cfg->ismomentum>0,
cfg->maxdetphoton*detreclen,cfg->seed,(uint)cfg->outputtype,0,0,cfg->faststep,
cfg->debuglevel,(uint)cfg->maxjumpdebug,cfg->gscatter,is2d};
if(param.isatomic)
Expand Down Expand Up @@ -1848,6 +1852,8 @@ void mcx_run_simulation(Config *cfg,GPUInfo *gpu){
#endif
if(cfg->issavedet)
sharedbuf+=gpu[gpuid].autoblock*sizeof(float)*(cfg->medianum-1);
if(cfg->ismomentum)
sharedbuf+=gpu[gpuid].autoblock*sizeof(float)*(cfg->medianum-1);
#ifdef USE_MT_RAND
sharedbuf+=(N+2)*sizeof(uint); // MT RNG uses N+2 uint in the shared memory
#endif
Expand Down
1 change: 1 addition & 0 deletions src/mcx_core.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ typedef struct __align__(16) KernelParams {
unsigned int issaveseed; /**< flag if one need to save the detected photon seeds for replay */
unsigned int issaveexit; /**< flag if one need to save the detected photon positions and dir vectors */
unsigned int issaveref; /**< flag if one need to save diffuse reflectance data in the 0-voxel layer next to the boundary */
unsigned int ismomentum; /**< 1 to save momentum transfer for detected photons, implies issavedet=1*/
unsigned int seedoffset; /**< offset of the seed, not used */
int seed; /**< RNG seed passted from the host */
unsigned int outputtype; /**< Type of output to be accummulated */
Expand Down
11 changes: 9 additions & 2 deletions src/mcx_utils.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@

const char shortopt[]={'h','i','f','n','t','T','s','a','g','b','B','z','u','H','P','N',
'd','r','S','p','e','U','R','l','L','-','I','-','G','M','A','E','v','D',
'k','q','Y','O','F','-','-','x','X','-','-','\0'};
'k','q','Y','O','F','-','-','x','X','-','-','m','\0'};

/**
* Long command line options
Expand All @@ -92,7 +92,8 @@ const char *fullopt[]={"--help","--interactive","--input","--photon",
"--printgpu","--root","--gpu","--dumpmask","--autopilot",
"--seed","--version","--debug","--voidtime","--saveseed",
"--replaydet","--outputtype","--outputformat","--maxjumpdebug",
"--maxvoidstep","--saveexit","--saveref","--gscatter","--mediabyte",""};
"--maxvoidstep","--saveexit","--saveref","--gscatter","--mediabyte",
"--momentum",""};

/**
* Output data types
Expand Down Expand Up @@ -202,6 +203,7 @@ void mcx_initcfg(Config *cfg){
cfg->debuglevel=0;
cfg->issaveseed=0;
cfg->issaveexit=0;
cfg->ismomentum=0;
cfg->replay.seed=NULL;
cfg->replay.weight=NULL;
cfg->replay.tof=NULL;
Expand Down Expand Up @@ -1675,6 +1677,10 @@ void mcx_parsecmd(int argc, char* argv[], Config *cfg){
case 'd':
i=mcx_readarg(argc,argv,i,&(cfg->issavedet),"char");
break;
case 'm':
i=mcx_readarg(argc,argv,i,&(cfg->ismomentum),"char");
if (cfg->ismomentum) cfg->issavedet=1;
break;
case 'r':
i=mcx_readarg(argc,argv,i,&(cfg->respin),"int");
break;
Expand Down Expand Up @@ -2011,6 +2017,7 @@ where possible parameters include (the first value in [*|*] is the default)\n\
-X [0|1] (--saveref) 1 to save diffuse reflectance at the air-voxels\n\
right outside of the domain; if non-zero voxels\n\
appear at the boundary, pad 0s before using -X\n\
-m [0|1] (--momentum) 1 to save photon momentum transfer,0 not to save\n\
-q [0|1] (--saveseed) 1 to save photon RNG seed for replay; 0 not save\n\
-M [0|1] (--dumpmask) 1 to dump detector volume masks; 0 do not save\n\
-H [1000000] (--maxdetphoton) max number of detected photons\n\
Expand Down
3 changes: 2 additions & 1 deletion src/mcx_utils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
@brief MCX configuration header
*******************************************************************************/

#ifndef _MCEXTREME_UTILITIES_H
#define _MCEXTREME_UTILITIES_H

Expand Down Expand Up @@ -167,6 +167,7 @@ typedef struct MCXConfig{
char issaveseed; /**<1 save the seed for a detected photon, 0 do not save*/
char issaveexit; /**<1 save the exit position and dir of a detected photon, 0 do not save*/
char issaveref; /**<1 save diffuse reflectance at the boundary voxels, 0 do not save*/
char ismomentum; /**<1 to save momentum transfer for detected photons, implies issavedet=1*/
char srctype; /**<0:pencil,1:isotropic,2:cone,3:gaussian,4:planar,5:pattern,\
6:fourier,7:arcsine,8:disk,9:fourierx,10:fourierx2d,11:zgaussian,12:line,13:slit*/
char outputtype; /**<'X' output is flux, 'F' output is fluence, 'E' energy deposit*/
Expand Down

0 comments on commit 63bce2b

Please sign in to comment.