MCXLAB Stalls when running high photon counts

[JCortes21799]

MCXLAB simulations hang at random intervals when running higher number of photons. Hanging has been seen from 1e7 to 1e10 photons. This has been observed with different sets of code and across machines. Machine that code has stalled the most has the following specifications:
Processor - AMD EPYC 7502P 32-Core
Linux Mint 20.1
CUDA driver - CUDA 11.2
GPU - GeForce GTX 1080

[fangq]

@JCortes21799, I will need your script/cfg in order to reproduce this first, we don't have Mint, but have Ubuntu 20.04 servers with cuda 11.2, 2080S

[JCortes21799]

This is the script/cfg being used. The CFG has also been run in individual cases and often works, however does hang on 1e9 and 1e10.

On Fri, Mar 5, 2021 at 11:43 AM Qianqian Fang ***@***.***> wrote: @JCortes21799 <https://github.com/JCortes21799>, I will need your script/cfg in order to reproduce this first, we don't have Mint, but have Ubuntu 20.04 servers with cuda 11.2, 2080S — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#114 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANAGD3AJKY4PVBXFBUM2UEDTCEC4PANCNFSM4YVQV4BQ> .

``` %Source Jentered Simulations %60 x 60 x 60 %Single channel %Changing the source to channel (s2c) distances for different trials %[depths], and diameters close all; clear all; seed1 = "623F9A9E"; seed2 = "fd0dd689"; seed3 = "12064fa6"; seed4 = "496edeaa"; diameters = [4 5 6]; NUM_TRIALS = 18; %depths 18 NUM_S2C = 26;%26 NUM_DIAM = 3; depths = [4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21]; %mm %diameter 2, diameter 1 done already, original case %need min_depth - diameter > 2? because layer 1 is air s2cs = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25] ; %% for dn = 1:NUM_DIAM for i=1:NUM_TRIALS for j=1:NUM_S2C % first query if you have supported GPU(s) channel_center = depths(i); dist_s2c = s2cs(j); diameter = diameters(dn); % define the simulation using a struct cfg.nphoton= 1e3; cfg.bc='cc-cc-'; cfg.seed=hex2dec(seed1); cfg.vol=uint8(ones(60,60,60)); z_start = (channel_center - diameter/2); z_end = (channel_center + diameter/2); channel_y = 30; cfg.vol( (channel_y-1) : (channel_y+1) , 1:60, z_start : z_end ) = 2; % add an inclusion cfg.prop=[0 0 1 1;0.005 1 0 1.37; 0.2 10 0.9 1.37]; % [mua,mus,g,n] cfg.issrcfrom0=1; src_y = dist_s2c + 29.5; cfg.srcpos=[src_y src_y 1]; cfg.srcdir=[0 0 1]; cfg.vol(:,:,1)=0; % pad a layer of 0s to get diffuse reflectance cfg.issaveref=1; cfg.gpuid=1; cfg.autopilot=1; cfg.tstart=0; cfg.tend=5e-9; cfg.tstep=5e-10; cfgs(1) = cfg; cfgs(2) = cfg; cfgs(1).prop = [0 0 1 1; 0.011 14.1 0.9 1.37 0.011 14.1 0.9 1.37]; cfgs(2).prop = [0 0 1 1; 0.011 14.1 0.9 1.37 0.21 14.1 0.9 1.37]; % calculate the fluence distribution with the given config [fluences,detpts,vols,seeds,trajs]=mcxlab(cfgs); % integrate time-axis (4th dimension) to get CW solutions cwdref_1=sum(fluences(1).dref,4); % diffuse reflectance cwdref_2=sum(fluences(2).dref,4); % diffuse reflectance % plot configuration and results %diff = squeeze(log(cwdref_1(:,:,1))) - squeeze(log(cwdref_2(:,:,1))); ch_img = cwdref_2(:,:,1); bg_img = cwdref_1(:,:,1); diff = ch_img - bg_img; diff = diff/bg_img; tr = num2str(channel_center); s1 = strcat('diffuse refle. at z=1; d = ', tr); num = num2str(src_y); s2 = strcat('; sp = ', num); title_string = strcat(s1, s2); % title(title_string); trial_str = num2str(channel_center) num_str = num2str(dist_s2c) diameter_str = num2str(diameter); seed_str = num2str(1); save_str_bg = strcat('seed_refined', '_', seed_str, '_', diameter_str, '_', trial_str, '_', num_str, '_', 'bg'); save_str_ch = strcat('seed_refined', '_', seed_str, '_',diameter_str, '_', trial_str, '_', num_str, '_', 'ch'); save_str_diff = strcat('seed_refined', '_', seed_str, '_',diameter_str, '_', trial_str, '_', num_str, '_', 'diff'); save(save_str_diff, 'diff') save(save_str_ch, 'ch_img') save(save_str_bg, 'bg_img') end end end ```

[fangq]

I ran a subset of your simulations on my ubuntu 20.04 box, it runs without any problem.

My suspicion is that your GPU is overheating after running a large number of simulation, if you don't have a good cooling system.

Next time, when you see this hanging again, please go to the terminal, and type nvidia-smi to see the temperature of the GPU(s). You can also install the lm-sensors package and type sensors command to see the CPU/PCI bus temperature.

also, it appears to me that you only need the fluence output, please remove all other unwanted outputs to speed up your computation. The output of detps and trajs are particularly memory/time consuming.

[fangq]

@JCortes21799, saw your offline message - from the nvidia-smi output, the GPU temperature looks fine.

I suspect that this is related to the progress bar printing issue we had previously

https://forums.developer.nvidia.com/t/using-pinned-memory-to-pass-simple-data-fails-on-windows/51351

and similarly for OpenCL

https://forums.developer.nvidia.com/t/solved-how-to-update-host-memory-variable-from-device-during-opencl-kernel-execution/59409/6

can you disable progress bar printing by setting cfg.debuglevel=''; or remove debuglevel from cfg, and rerun your simulation?

[JCortes21799]

Good evening, I reran the simulations without the progress bar and I still am seeing the same stalling issue. The nvidia-smi output was normal as well. This was the output across two different instances of AWS p3.2xlarge machines. Best, Juan

…

On Tue, Aug 24, 2021 at 3:08 PM Qianqian Fang ***@***.***> wrote: @JCortes21799 <https://github.com/JCortes21799>, saw your offline message - from the nvidia-smi output, the GPU temperature looks fine. I suspect that this is related to the progress bar printing issue we had previously https://forums.developer.nvidia.com/t/using-pinned-memory-to-pass-simple-data-fails-on-windows/51351 and similarly for OpenCL https://forums.developer.nvidia.com/t/solved-how-to-update-host-memory-variable-from-device-during-opencl-kernel-execution/59409/6 can you disable progress bar printing by setting cfg.debuglevel=''; or remove debuglevel from cfg, and rerun your simulation? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#114 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANAGD3E4YUXT4CGFQ567MQ3T6PUZBANCNFSM4YVQV4BQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .

[fangq]

@JCortes21799, the sample script you provided above runs 3x18x26 simulations, does the stalling always happen at a specific i/j/dn combination or random?

if it happens at a specific combination, then please set a break-point for that failed mcxlab call and save the cfg input into a .mat file and send it to me, I need to first reproduce this in order to debug.

if it happens randomly, I strongly suspect that this is related to your hardware/network stability, please consult your hardware provider to debug the issue.

[JCortes21799]

The code seems to hang always on j = 1 so I have attached the CFG's from code ran on four different AWS EC2 p3.2xlarge instances as well. I have also added the updated script I have been using.
bugged_cfgs.zip

This is the updated script I am running:

close all;
clear all;

seed1   = "623F9A9E";

SCALE = 1; %10x per side
NUM_TRIALS = 1; %depths 18 4
NUM_S2C    = 2; %source-to-channel distances 26 10
NUM_DIAM   = 1;  %diameters 3

X_DIM = 31; %60 * SCALE + 1 for the centering
Y_DIM = 61; %60 * SCALE + 1 for the centering
Z_DIM = 31; %60 * SCALE + 1 for the centering
CENTER_Y = (Y_DIM - 1)/2 + 1; %301 -> Value is 31
CENTER_X = (X_DIM - 1)/2 + 1; %151 -> Value is 16

depths = [5 10 15 20] * SCALE; %mm
%diameter 2, diameter 1 done already, original case
%need min_depth - diameter > 2? because layer 1 is air
s2cs = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25] * SCALE;
 %%
    for i=1:NUM_TRIALS
        for j=1:NUM_S2C
            channel_center = depths(i); %SCALED
            dist_s2c = s2cs(j); %SCALED
            diameter = 4*SCALE; %SCALED
            % define the simulation using a struct
            cfg.nphoton= 1e9;
            cfg.bc='cc-cc-';
            cfg.seed=hex2dec(seed1);
            cfg.vol=ones(Y_DIM, X_DIM, Z_DIM); %SCALED 

            channel_y = CENTER_Y; %This is now perfectly centered

            %Build a slice that has the cylinder circle as 2's and the rest as 1's
            slice = ones(Z_DIM, Y_DIM); %Row is Z, Col is Y for the phantom
            radius = diameter/2;
            for k = 1:Y_DIM %y
                for L = 1:Z_DIM  %z
                    x1 = CENTER_Y;
                    y1 = channel_center;
                    x2 = k;
                    y2 = L;
                    distance = sqrt( (x1 - x2)^2 + (y1 - y2)^2 );
        
                    if(distance <= radius)
                        slice(L, k) = 2;
                    end
                end
            end
        
            %Compound many slices into the phantom volume
            for slice_i =1:X_DIM
                cfg.vol(:, slice_i, :) = slice';
            end
            
            %Need to pick the tissue PROPERTIES *****************

            %Air, Surrounding Tissues, Inclusion = [0 ;1 ;2] (inclusion is the vessel filled with blood)
            cfg.prop=[0 0 1 1;
            0.022 0.47 0.9 1.4;
            0.022 0.47 0.9 1.4]; % [mua,mus,g,n]
            cfg.issrcfrom0=1;             %1-first voxel is [0 0 0], [0]- first voxel is [1 1 1]
            src_y = dist_s2c + channel_y; %Moving Y for the Source
            src_x = CENTER_X;               %301 %Stationary X for the Source
            cfg.srcpos=[src_y src_x 1];   %Only Y is moving
            cfg.srcdir=[0 0 1];           %Downwards
            cfg.vol(:,:,1)=0;             %Pad a layer of 0s to get diffuse reflectance, look at properties above
            cfg.issaveref=1;
            cfg.gpuid=1;
            cfg.autopilot=1;
            cfg.tstart=0;
            cfg.tend=5e-9;
            cfg.tstep=5e-9;
	    cfg.unitinmm = 0.1;
        %cfg.respin = 2;
        %cfg.debuglevel = 'P'; 
	    %cfg.workload = [60]; 
        %Make two cfg copies that are the same
        cfgs(1) = cfg;
        cfgs(2) = cfg;

        %First is case where channel has same properties as surrounding tissue [no blood flowing]
        cfgs(1).prop = [0 0 1 1;
        0.02 12 0.9 1.4;
        0.02 12 0.9 1.4];

        %Second is case where channel has higher absorption properties [blood flowing]
        cfgs(2).prop = [0 0 1 1;
        0.02  12 0.9 1.4;
        0.2   17 0.9 1.4];

        disp("----------------------Configurations Start Here----------------------");
        cfgs
        disp("----------------------Configuration Ends Here----------------------");
    
        disp("----------------------Configuration 1 Start Here----------------------");
        cfgs(1)
        disp("----------------------Configuration 1 Ends Here----------------------");
            
        disp("----------------------Configuration 2 Start Here----------------------");
        cfgs(2)
        disp("----------------------Configuration 2 Ends Here----------------------");
        
        %Save the configurations
        bug_base_str = 'bugged_cfgs';
        i_str = num2str(i);
        j_str = num2str(j);
        bug_save_str = strcat(bug_base_str, '_', i_str, '_', j_str);
        save(bug_save_str, 'cfgs');
        
        
        % calculate the fluence distribution with the given config
        fluences=mcxlab(cfgs);

        % integrate time-axis (4th dimension) to get CW solutions
        cwdref_1=sum(fluences(1).dref,4);     % diffuse reflectance
        cwdref_2=sum(fluences(2).dref,4);     % diffuse reflectance

        %Take the top slices of the resulting data, what the Camera sees
        ch_img = cwdref_2(:,:,1); 
        bg_img = cwdref_1(:,:,1);
        diff = ch_img - bg_img;
        diff = diff/bg_img;

        %Strings for plots
        tr = num2str(channel_center);
        s1 = strcat('diffuse refle. at z=1; d = ', tr);
        num = num2str(src_y);
        s2 = strcat('; sp = ', num);
        title_string = strcat(s1, s2);
%           title(title_string);

        %Saving the images for the trials
        trial_str = num2str(channel_center/SCALE)
        num_str = num2str(dist_s2c/SCALE)
        diameter_str = num2str(diameter/SCALE)
        seed_str = num2str(1);
        base_string = 'aws_small_1e9_progress'

        save_str_bg   = strcat(base_string, '_', seed_str, '_', diameter_str, '_', trial_str, '_', num_str, '_', 'bg');
        save_str_ch   = strcat(base_string, '_', seed_str, '_', diameter_str, '_', trial_str, '_', num_str, '_', 'ch');
        save_str_diff = strcat(base_string, '_', seed_str, '_', diameter_str, '_', trial_str, '_', num_str, '_', 'diff');
        save(save_str_diff, 'diff')
        save(save_str_ch, 'ch_img')
        save(save_str_bg, 'bg_img')
        end
    end
 

[fangq]

@JCortes21799, I was able to reproduce the hanging from your above code. The issue appears to come from the use of cyclic bc in your line cfg.bc='cc-cc-'. if you disable that line, the script runs fine.

I will try to find some time to look into this. in the meantime, I suggest you disable this line because it may not do what you wanted to achieve - if cyclic bc is used, essentially your source is mirrored in x/y directions to form an infinite x/y grid pencil beam array. So, it may not necessarily give you a better solution for simulating an infinite slab in x/y (where you just want to have a single source). If you make your cube wide enough, it should give you a rather good approximation without using this bc.