Skip to content

Commit

Permalink
[bug] enable photon sharing for pattern3d source
Browse files Browse the repository at this point in the history
  • Loading branch information
@fangq
fangq committed Nov 21, 2024
1 parent bd62362 commit dcf11c4
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 8 deletions.
21 changes: 15 additions & 6 deletions src/mcx_core.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1282,10 +1282,10 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*

p->w = 1.f;
}
} else if (gcfg->srctype == MCX_SRC_FOURIER)
} else if (gcfg->srctype == MCX_SRC_FOURIER) {
p->w = launchsrc->pos.w * (cosf((floorf(launchsrc->param1.w) * rx + floorf(launchsrc->param2.w) * ry
+ launchsrc->param1.w - floorf(launchsrc->param1.w)) * TWO_PI) * (1.f - launchsrc->param2.w + floorf(launchsrc->param2.w)) + 1.f) * 0.5f; //between 0 and 1
else if (gcfg->srctype == MCX_SRC_PENCILARRAY) {
} else if (gcfg->srctype == MCX_SRC_PENCILARRAY) {
p->x = launchsrc->pos.x + floorf(rx * launchsrc->param1.w) * launchsrc->param1.x / (launchsrc->param1.w - 1.f) + floorf(ry * launchsrc->param2.w) * launchsrc->param2.x / (launchsrc->param2.w - 1.f);
p->y = launchsrc->pos.y + floorf(rx * launchsrc->param1.w) * launchsrc->param1.y / (launchsrc->param1.w - 1.f) + floorf(ry * launchsrc->param2.w) * launchsrc->param2.y / (launchsrc->param2.w - 1.f);
p->z = launchsrc->pos.z + floorf(rx * launchsrc->param1.w) * launchsrc->param1.z / (launchsrc->param1.w - 1.f) + floorf(ry * launchsrc->param2.w) * launchsrc->param2.z / (launchsrc->param2.w - 1.f);
Expand Down Expand Up @@ -3223,12 +3223,13 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {

CUDA_ASSERT(cudaMemcpy(genergy, energy, sizeof(float) * (gpu[gpuid].autothread << 1), cudaMemcpyHostToDevice));

if (cfg->srcpattern)
if (cfg->srcpattern) {
if (cfg->srctype == MCX_SRC_PATTERN) {
CUDA_ASSERT(cudaMemcpy(gsrcpattern, cfg->srcpattern, sizeof(float) * (int)(cfg->srcparam1.w * cfg->srcparam2.w * cfg->srcnum), cudaMemcpyHostToDevice));
} else if (cfg->srctype == MCX_SRC_PATTERN3D) {
CUDA_ASSERT(cudaMemcpy(gsrcpattern, cfg->srcpattern, sizeof(float) * (int)(cfg->srcparam1.x * cfg->srcparam1.y * cfg->srcparam1.z * cfg->srcnum), cudaMemcpyHostToDevice));
}
}

/**
* Copy constants to the constant memory on the GPU
Expand Down Expand Up @@ -3695,12 +3696,16 @@ is more than what your have specified (%d), please use the -H option to specify
*/
#pragma omp master
{
if (cfg->issave2pt && cfg->srctype == MCX_SRC_PATTERN && cfg->srcnum > 1) { // post-processing only for multi-srcpattern
if (cfg->issave2pt && (cfg->srctype == MCX_SRC_PATTERN || cfg->srctype == MCX_SRC_PATTERN3D) && cfg->srcnum > 1) { // post-processing only for multi-srcpattern
srcpw = (float*)calloc(cfg->srcnum, sizeof(float));
energytot = (float*)calloc(cfg->srcnum, sizeof(float));
energyabs = (float*)calloc(cfg->srcnum, sizeof(float));
int psize = (int)cfg->srcparam1.w * (int)cfg->srcparam2.w;

if (cfg->srctype == MCX_SRC_PATTERN3D) {
psize = (int)cfg->srcparam1.x * (int)cfg->srcparam1.y * (int)cfg->srcparam1.z;
}

for (i = 0; i < int(cfg->srcnum); i++) {
float kahanc = 0.f;

Expand Down Expand Up @@ -3804,10 +3809,14 @@ is more than what your have specified (%d), please use the -H option to specify
/**
* In photon sharing mode, where multiple pattern sources are simulated, each solution is normalized separately
*/
if (cfg->srctype == MCX_SRC_PATTERN && cfg->srcnum > 1) { // post-processing only for multi-srcpattern
if ((cfg->srctype == MCX_SRC_PATTERN || cfg->srctype == MCX_SRC_PATTERN3D) && cfg->srcnum > 1) { // post-processing only for multi-srcpattern
float scaleref = scale[0];
int psize = (int)cfg->srcparam1.w * (int)cfg->srcparam2.w;

if (cfg->srctype == MCX_SRC_PATTERN3D) {
psize = (int)cfg->srcparam1.x * (int)cfg->srcparam1.y * (int)cfg->srcparam1.z;
}

for (i = 0; i < int(cfg->srcnum); i++) {
scale[i] = psize / srcpw[i] * scaleref;
}
Expand Down Expand Up @@ -3920,7 +3929,7 @@ is more than what your have specified (%d), please use the -H option to specify
((cfg->issavedet == FILL_MAXDETPHOTON) ? cfg->energytot : ((double)cfg->nphoton * ((cfg->respin > 1) ? (cfg->respin) : 1))) / max(1, cfg->runtime));
fflush(cfg->flog);

if (cfg->issave2pt && cfg->srctype == MCX_SRC_PATTERN && cfg->srcnum > 1) {
if (cfg->issave2pt && (cfg->srctype == MCX_SRC_PATTERN || cfg->srctype == MCX_SRC_PATTERN3D) && cfg->srcnum > 1) {
for (i = 0; i < (int)cfg->srcnum; i++) {
MCX_FPRINTF(cfg->flog, "source #%d total simulated energy: %.2f\tabsorbed: " S_BOLD "" S_BLUE "%5.5f%%" S_RESET"\n(loss due to initial specular reflection is excluded in the total)\n",
i + 1, energytot[i], energyabs[i] / energytot[i] * 100.f);
Expand Down
2 changes: 1 addition & 1 deletion src/mcx_utils.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2852,7 +2852,7 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
MCX_ERROR(-1, "Optode.Source.Pattern JData-annotated array must be in the 'single' format");
}

if (ndim == 3 && dims[2] > 1 && dims[0] > 1 && cfg->srctype == MCX_SRC_PATTERN) {
if (ndim == 3 && dims[2] > 1 && dims[0] > 1 && cfg->srctype == MCX_SRC_PATTERN || cfg->srctype == MCX_SRC_PATTERN3D) {
cfg->srcnum = dims[0];
}
} else {
Expand Down
6 changes: 5 additions & 1 deletion src/mcxlab.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1130,10 +1130,14 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
arraydim = mxGetDimensions(item);
dimtype dimz = 1;

if (mxGetNumberOfDimensions(item) == 3) {
if (mxGetNumberOfDimensions(item) >= 3) {
dimz = arraydim[2];
}

if (mxGetNumberOfDimensions(item) == 4) {
dimz *= arraydim[3];
}

double* val = mxGetPr(item);

if (cfg->srcpattern) {
Expand Down

0 comments on commit dcf11c4

Please sign in to comment.