Total reflection energy

[zhang haohui]

Dear Prof. Fang,

I am studying a simple problem using mcxlab: a pencil ray is injected into the center of a 15cm*15cm*7.5cm block. I tried to compute the total energy absorbed and reflected by the block. I expected their sum is 1. 

However, I found the absorbed energy = sum(CWfluence,'all')*0.01 = 0.18, and it is the same as printed text by the software. The total reflection is sum(cwdref(:,:,1),'all') is much smaller than 1-0.18. I have verified the energy escaped from all other five surfaces is zero. So the reflected energy I got is incorrect. I found that cwdref(150,150,1) is the same as its neighborhood pixels. However, it should be the largest because it is exactly beneath the light source. I guess the error is due to inaccurate reflected energy at this pixel (150,150,1).

I have attached my code below. Could you please help me check out the reason? Thank you very much!

% Bulid the shape

Vfinal = zeros(302,302,152);

Vfinal(2:301,2:301,2:151) = 1;

% Vfinal(:,:,end) = 0;

%% prepare cfg for MCX simulation

clear cfg

cfg.nphoton=1e9;

cfg.outputtype='fluence';

cfg.seed = 77542;

% tissue labels:0-ambient air,1-tissue

cfg.vol=Vfinal;

cfg.prop = [0,0,1,1;0.01,2,0,1];

% light source

cfg.srctype = 'pencil';

cfg.srcpos=[149.5,149.5,0];

cfg.srcdir=[0,0,1];

cfg.issrcfrom0=1;

% time windows

cfg.tstart=0;

cfg.tend=1e-8;

cfg.tstep=1e-9;

% other simulation parameters

cfg.isspecular=1;

cfg.isreflect=1;

cfg.autopilot=1;

cfg.gpuid=1;

cfg.issaveref=1; % save diffuse reflectance

cfg.unitinmm = 1;

% Detector

cfg.detpos=[50,50,100,2];

tic

%% run MCX simulation

[flux,detps]=mcxlab(cfg);

toc

%% post-simulation data processing and visualization

% convert time-resolved fluence to CW fluence

CWfluence=sum(flux.data,4);

cwdref=sum(flux.dref,4); % diffuse reflectance

absorb = sum(CWfluence,'all')*0.01;

reflect = sum(cwdref(:,:,1),'all');

[Qianqian Fang]

the two quantities you computed, absorb and reflect, have different units, one is fluence (J/mm^2), the other is energy (J).

you should change the below line

cfg.outputtype='fluence';

to

cfg.outputtype='energy';

after removing the *0.01, this should give you a unity output for absorb+reflect;

by the way, you don't need 1e9 photons to show this, even a small number of photon should be able to show this energy conservation.

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/9f6f43f4-4296-48d5-8e4d-4c22c74a7479n%40googlegroups.com.

[Qianqian Fang]

oh, forgot to mention, another reason that you are not getting conserved energy is because you only summed the dref along z=0 plane; photons escaped from other directions are not included.

you should change the summation to

absorb = sum(CWfluence(:))

reflect = sum(cwdref(:))

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/16cbbfd3-9b61-4565-360e-90847b50b30b%40neu.edu.

[zhang haohui]

Hello, Prof. Fang,

Thanks for your quick response. I've changed the outputtype according to your suggestion, however, the absorb+reflect is still not equal to 1. In fact, when cfg.nphoton=1e6 or 1e7, the absorb+reflect =1. When cfg.nphoton=1e6 or 1e7, the absorb+reflect<1. By comparing the cwdref of 1e6 and 1e8, I found cwdref(150,150,1) has a significant difference (0.31 v.s. 0.17) and led to the overall difference in absorb+reflect. The summation of reflect -  cwdref(150,150,1)

Here I wrote my revised code.

% Bulid the shape

Vfinal = zeros(302,302,152);

Vfinal(2:301,2:301,2:151) = 1;

% Vfinal(:,:,end) = 0;

%% prepare cfg for MCX simulation

clear cfg

cfg.nphoton=1e6;

cfg.outputtype='energy';

absorb = sum(CWfluence,'all');

reflect = sum(cwdref(:,:,:),'all');

absorb + reflect - cwdref(150,150,1)

absorb + reflect

Thanks for your help.

Best wishes

Haohui

[Qianqian Fang]

Haohui,

thanks, now I see the issue. I also understood now why you designed such a large domain - as you intended to minimize the escaping photons from other sides.

I was able to reproduce the loss of total energy (absorb+reflect) in 1e8 and 1e9 cases. I recognize that this is a bug in the code - reminding me that I had encountered a similar bug previously.

I just created a ticket for this problem and explained the reasons, please see #195

https://github.com/fangq/mcx/issues/195

basically, this is related to another high-priority bug I fixed back in 2018 (#41)

https://github.com/fangq/mcx/issues/41

what happened was that mcx uses "single-precision" floating point numbers for all fluence calculations instead of "double precision". This is because modern GPUs primarily uses single-precision float (fp32) for the bulk of the computation; double-precision hardware is very weak (1-vs-32 compared to single hardware) on consumer-grade GPUs (gaming or laptop GPUs), and only has strong presence on enterprise grade GPUs (Tesla), which are very expensive.

Because of this, we have optimized mcx to avoid double-precision (fp64) math or storage and exclusively using fp32. However, fp32 has limitations when performing numerous accumulations like the one used in mcx: when adding a small energy loss over a large existing number, such as a voxel near the source, the round-off error is quite large and resulting in loss of precision (and effectively loss of energy).

In #41, we developed a double buffer approach to avoid such loss using pure fp32 storage (see the figures in #41), however, it was only applied to fluence, but not applied to dref accumulation.

the issue you reported is related to such loss in dref (i.e. your reflect value) due to the above reason. I just applied the following patch to fix this

https://github.com/fangq/mcx/commit/9220578840024b296843dc595432471086037a25

after using this patch, the loss in 1e8 an 1e9 are gone. the absorb+reflect values are reasonably close to 1.

I just recompiled all mcx-related nightly-build packages, please download from https://mcx.space/nightly/ and test again. I have also updated pmcx to v0.2.6 to include the fix to this critical bug.

by the way: if you know how to recompile mcx, you can run "make double" or "make moredouble" to use double-precision buffer to accumulate fluence/dref. This is an alternative solution, but its speed depends on your GPU fp64 hardware.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/9d805f87-9ff1-4bc4-bd37-45802e18dc04n%40googlegroups.com.

[zhang haohui]

Hello, Prof. Fang,

Thank you very much! I just tested your new mcxlab code and the total energy is 1 for 1e8 and 1e9 photons! It works great!

Best wishes

Haohui