disordered-photonics / celes Goto Github PK

Allow the user to switch between single and double precision.

Dear authors,

This code seems to be not able to solve "periodic particle array" problem with "plane wave sources".

Do you have any updated version of Celes to add a "periodic particle array" boundary condition inside code?

Thank you!

Le

The try-catch logic in compute_maximal_particle_distance (celes_particles.m) doesn't make sense. Got the error: "Undefined function or variable 'CH'."

While trying to understand the framework for my Master-Thesis, I noticed a frequent use of for-loops throughout the code.
Usually such loops are quite inefficient and should be swaped out for matrix/vector multiplications.
Is there and intent behind the loops or was it used to just solve the problem?
I am asking since I wanted to try to improve the performance by removing such for-loops, which can be swaped out, but I wanted to make sure I am not doing it in vain.

A method that can enhance the efficiency of CELES a lot is described in this paper by Chew.

Implementing this would be a bigger project, but should be doable.

Animations (GIFS) always impress. We could add export functionality to create movies from the field plots.

I have modified the examples and change the distribution of spheres. I am prompted to use reshape incorrectly.
And my txt of spheres was shown below.
1.txt

And this is my main code.
`% =========================================================================
%> @brief Run this script to start the simulation
% =========================================================================

%% add all folders to the MATLAB search path
addpath(genpath('./src'))

%% import example file with sphere positions, radii and refractive indices
data = dlmread('examples/1.txt');

%% initialize the CELES class instances

% initialize particle class instance
% - positionArray: Nx3 array (float) in [x,y,z] format
% - refractiveIndexArray: Nx1 array (complex) of complex refractive indices
% - radiusArray: Nx1 array (float) of sphere radii
particles = celes_particles('positionArray', data(:,1:3), ...
'refractiveIndexArray', data(:,5)+1i*data(:,6), ...
'radiusArray', data(:,4) ...
);

% initialize initial field class instance
% - polarAngle: scalar (float) polar angle of incoming beam/wave,
% in radians. for Gaussian beams, only 0 or pi are
% currently possible
% - azimuthalAngle: scalar (float) azimuthal angle of incoming
% beam/wave, in radians
% - polarization: string (char) polarization of incoming beam/wave
% ('TE' or 'TM')
% - beamWidth: scalar (float) width of beam waist. use 0 or inf
% for plane wave
% - focalPoint: 1x3 array (float) focal point
initialField = celes_initialField('polarAngle', 0, ...
'azimuthalAngle', 0, ...
'polarization', 'TE', ...
'beamWidth', 2000, ...
'focalPoint', [0,0,0] ...
);

% initialize input class instance
% - wavelength: scalar (float) vacuum wavelength, same unit as
% particle positions and radii
% - mediumRefractiveIndex: scalar (complex) refractive index of environment
% - particles: valid instance of celes_particles class
% - initialField: valid instance of celes_initialField class
input = celes_input('wavelength', 550, ...
'mediumRefractiveIndex', 1.49, ...
'particles', particles, ...
'initialField', initialField ...
);

% initialize preconditioner class instance
% - type: string (char) type of preconditioner (currently
% only 'blockdiagonal' and 'none' available)
% - partitionEdgeSizes 1x3 array (float) edge size of partitioning cuboids
% (applies to 'blockdiagonal' type only)
precnd = celes_preconditioner('type', 'blockdiagonal', ...
'partitionEdgeSizes', [1200,1200,1200] ...
);

% initialize solver class instance
% - type: string (char) solver type (currently 'BiCGStab' or
% 'GMRES' are implemented)
% - tolerance: scalar (float) target relative accuracy of solution
% - maxIter: scalar (int) maximum number of iterations allowed
% - restart: scalar (int) restart parameter (applies only to
% GMRES solver)
% - preconditioner: valid instance of celes_preconditioner class
solver = celes_solver('type', 'GMRES', ...
'tolerance', 5e-4, ...
'maxIter', 1000, ...
'restart', 200, ...
'preconditioner', precnd);

% initialize numerics class instance
% - lmax: scalar (int) maximal expansion order of scattered
% fields around particle center
% - polarAnglesArray: 1xN array (float) sampling of polar angles in the
% plane wave patterns, in radians
% - azimuthalAnglesArray: sampling of azimuthal angles in the plane wave
% patterns, in radians
% - gpuFlag: scalar (bool) set to false if you experience GPU
% memory problems at evaluation time (translation
% operator always runs on GPU, though)
% - particleDistanceResolution: scalar (float) resolution of lookup table for
% spherical Hankel function (same unit as wavelength)
% - solver: valid instance of celes_solver class
numerics = celes_numerics('lmax', 3, ...
'polarAnglesArray', 0:pi/5e3:pi, ...
'azimuthalAnglesArray', 0:pi/1e2:2*pi, ...
'gpuFlag', true, ...
'particleDistanceResolution', 1, ...
'solver', solver);

% define a grid of points where the field will be evaulated
[x,z] = meshgrid(-5000:50:5000, -5000:50:5000); y = zeros(size(x));
% initialize output class instance
% - fieldPoints: Nx3 array (float) points where to evaluate the
% electric near field
% - fieldPointsArrayDims: 1x2 array (int) dimensions of the array, needed to
% recompose the computed field as a n-by-m image
output = celes_output('fieldPoints', [x(:),y(:),z(:)], ...
'fieldPointsArrayDims', size(x));

% initialize simulation class instance
% - input: valid instance of celes_input class
% - numerics: valid instance of celes_input class
% - output: valid instance of celes_output class
simul = celes_simulation('input', input, ...
'numerics', numerics, ...
'output', output);

%% run simulation
simul.run;

% evaluate transmitted and reflected power
simul.evaluatePower;
fprintf('transmitted power: \t%.4f %%\n', ...
simul.output.totalFieldForwardPower/simul.output.initialFieldPower100)
fprintf('reflected power: \t%.4f %%\n', ...
simul.output.totalFieldBackwardPower/simul.output.initialFieldPower100)

% evaluate field at output.fieldPoints
simul.evaluateFields;

%% plot results
% display particles
figure('Name','Particle positions','NumberTitle','off');
plot_spheres(gca,simul.input.particles.positionArray, ...
simul.input.particles.radiusArray, ...
simul.input.particles.refractiveIndexArray)

% plot near field
figure('Name','Near-field cross-cut','NumberTitle','off');
plot_field(gca,simul,'abs E','Total field')
colorbar
caxis([0,2])

% % export animated gif
% figure('Name','Animated near-field cross-cut','NumberTitle','off');
% plot_field(gca,simul,'real Ey','Total field','Ey_total.gif')
`

Dear colleagues,
I am trying to run your exciting code. I got the following message when cloning celes as prescribed in readme.md:

git clone [email protected]:disordered-photonics/celes.git
Cloning into 'celes'...
Permission denied (publickey).
fatal: Could not read from remote repository.

Please make sure you have the correct access rights
and the repository exists.

Did I miss a step?
Best regards,

Guillaume Legros
http://www.dalembert.upmc.fr/home/legros/

The code has a margin to improve the performance, but I saw that the PR and issue response had been stopped since 2020.
Would you be able to take my potential PRs and update this solver?

right now, the user accesses the polarization property ("TE" or "TM") and the hidden pol property is automatically set during intialization. The problem is, that if the polarization is changed later, one needs to call the setPolIdx method again. It would be good if this happens automatically.

I would suggest to change the pol property to a dependent property and rename the setPolIdx to get.pol method. Then, pol is updated when polarization is updated. @lpattelli, does this conflict with how you intended the initialization to happen?

I'd like to suggest the addition of the capability to run models with multiple spheres arbitrarily placed inside other spheres (similar to MSTM). As far as I can tell, there is not currently another light scattering code with this capability that uses GPU computing, and it would represent a significant increase in capability.

I'm happy to donate my time and coding skills, but unfortunately don't currently know enough to do it myself.

Running the example file CELES_MAIN.m with a non-zero imaginary part of 'mediumRefractiveIndex' causes the execution to terminate with an error right after it says
"compute scattered field plane wave pattern:"
The error is:
"Error using celes_planeWavePattern/validatePropertiesImpl (line 123) invalid k value: Expected input to be real."
the refractive index I used was 1.2 + 0.01i and the light wavelength 500. Also, output.initialFieldPower is complex, with large imaginary part, in my case 3.7e6 + 3.1e4i.
Does anyone else have this issue?

Hello, I am trying to decide which nvidia gpu to buy.
Roughly how many particles can I simulate using 3GB gpu ram with default settings ?

In the function computeTotalFieldPWP of the celes_simulation class, obj.output.scatteredFieldPlaneWavePattern is overwritten by obj.output.totalFieldPlaneWavePattern. Maybe it is another consequence of celes_output deriving from Matlab.System.

Right now, the user needs to have a setup with Matlab, CUDA and C++ compiler available to even test our software.

It would be great to have an alternative Matrix-Vector product that runs on the CPU only. This function will probably be slower than the current, but at least allow the user to test the software.

This could be achieved either by a Matlab script coupling_matrix_multiply_CPU.m or by a c++ file that does what coupling_matrix_multiply_CUDA.cu does (and which is compiled by mex).

Up to now, CELES simulations were saved by trivially saving the celes_simulation object as a whole. This seems broken since recent updates in the data structure (Lorenzo, is that correct?).
Another drawback of saving the whole celes_simulation object is that lookup tables and other heavy stuff might be saved leading to large file sizes.

To meet this issue, we could implement a save method in the simulation class that manually stores the relevant attributes into some data structure which is then saved, as well as a load method that reads from that structure and initializes a simulation object.

Dear Drs. Egel and Pattelli,

I learned about CELES through your 2017 publication titled ‘ CELES: CUDA-accelerated simulation of electromagnetic scattering by large ensembles of spheres.’ Since my research involves studying the diffraction patterns generated by organized particles, I am very interested in using CELES to confirm some of my experimental outcomes.

Over the past week, I was able to install CELES correctly and ran some trial simulations. However, I do have a few questions about CELES, which I hope you can answer for me.

(1) Can CELES be used for a single scatterer? I want to perform a comparison between CELES and Mie Scattering theory.
(2) Is the unit of all dimensional parameters (e.g., particle radius) nm?
(3) Is it possible to implement a collimated beam (not focused) with a finite diameter? If so, could you please advise how I may do it?
(4) Is there an up limit regarding the number of scatterers? In your article (figure 6), you mentioned it is possible to simulate 10^5 particles.
(5) Your article mentioned that CELES could simulate the E field in the far-field. Would it be possible to study the scattered E field, for example, 10 cm away from the scatterer cluster?

I look forward to hearing from you soon.

WC Lin
[email protected]
Florida International University

I am currently analyzing the behavior of particle clusters under different wavelengths. For this, I need the expansion coefficients a and f, which are provided by the tables properties of the simulation class (initialFieldCoefficients and scatteredFieldCoefficients). Since celes can only accept a single wavelength, I need to rerun the same code multiple times, which does take a lot of time. I noticed that the run function of the simulation class accepts an initial guess of the scattered field coefficients f0, so I am using the coefficients of the previous run to initialize them. My reason behind this is that the field should not differ much by changing the wavelength slightly. I have not found similar parameters where I could provide some initial guesses in the preconditioner class, which takes also quite a long time to finish.
I was wondering if someone could advise me on how I could speed up things.
Thanks for the help!

The method compute_maximal_particle_distance calls convhull which raises an error if the particle number is below 4.

Hi!
Thanks for developing this! Due to the SSD failure, our desktop that was able to run this program crashed, and we fail to retrive the relevant information related to environment and software version. We tried to rebuild the environment, but in newer version of matlab and cuda (2021a with CUDA 10.2, 11), we get this error message

number of particles: 500
compiling CUDA-code with lmax=3.
Building with 'NVIDIA CUDA Compiler'.
Error using mex
ptxas fatal   : Unresolved extern function '_Z3powfi'
nvcc warning : The 'compute_35', 'compute_37', 'compute_50', 'sm_35', 'sm_37' and 'sm_50' architectures are deprecated, and may be removed in a future release
(Use -Wno-deprecated-gpu-targets to suppress warning).
coupling_matrix_multiply_CUDA.cu



Error in mexcuda (line 168)
    [varargout{1:nargout}] = mex(mexArguments{:});

Error in cuda_compile (line 73)
    mexcuda(verb_flag, DEFINES, '-outdir', script_dir, '-output', OUT_FILE, SRC_FILE)

Error in celes_simulation/run (line 394)
                cuda_compile(obj.numerics.lmax);

Error in CELES_MAIN (line 114)
simul.run;

We also tried the combination of 2017b, cuda8.0 and visual studio 15, but some error related to compiler pops up

Error using mex
No supported compiler or SDK was found. For options, visit https://www.mathworks.com/support/compilers.
Error in mexcuda (line 157)
    [varargout{1:nargout}] = mex(mexArguments{:});
Error in cuda_compile (line 73)
    mexcuda(verb_flag, DEFINES, '-outdir', script_dir, '-output', OUT_FILE, SRC_FILE)
Error in celes_simulation/run (line 394)
                cuda_compile(obj.numerics.lmax);
Error in CELES_MAIN (line 114)
simul.run;

Does anyone have ideas how to solve it? Thanks!

The SMARTIES software allows to compute the T-matrix of spheroids, which could be input to CELES. Dominik has already worked out how to extract the T-matrix from SMARTIES and convert it to the normalization that is used in CELES, he could probably forward his matlab files.

In Alan's forked version of this project, there also exists a branch that supports T-matrices other than spheres (from NFM-DS in his case).

The GUI is broken since d0263ba. In order to fix it, the routines celes/src/gui/app2simulation.m
and celes/src/gui/simulation2app.m must be adapted to the new constructor and data structure of the celes_simulation class.