From dcf11c4dfb8fdb38f830b4372c404f2c7541f5e5 Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Thu, 21 Nov 2024 00:41:24 -0500 Subject: [PATCH] [bug] enable photon sharing for pattern3d source --- src/mcx_core.cu | 21 +++++++++++++++------ src/mcx_utils.c | 2 +- src/mcxlab.cpp | 6 +++++- 3 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 75056ecd..9709c4b5 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -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); @@ -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 @@ -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; @@ -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; } @@ -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); diff --git a/src/mcx_utils.c b/src/mcx_utils.c index 74818b4f..699612f9 100644 --- a/src/mcx_utils.c +++ b/src/mcx_utils.c @@ -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 { diff --git a/src/mcxlab.cpp b/src/mcxlab.cpp index 38d146e4..df58c211 100644 --- a/src/mcxlab.cpp +++ b/src/mcxlab.cpp @@ -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) {