docs updates

docs/Options/performance.rst

@@ -0,0 +1,41 @@
+Performance
+===========
+
+Ray-tracing and RCWA are both computationally intensive. RayFlare currently has two approaches for
+increasing the speed of calculations:
+
+- **Parallel computation**: RayFlare uses the joblib library to parallelize computations (over CPUs) where possible.
+  In ray-tracing, the parallelisation is done over rays (one wavelength calculated at a time, multiple rays - note that
+  prior to v2.0.0 the parallelisation was over wavelengths instead). In RCWA, the parallelisation is over wavelengths.
+- **Numba JIT (just-in-time) compilation**: RayFlare uses the Numba library to compile many ray-tracing functions to machine
+  code, which can result in significant speed-ups.
+
+By default, both of these options are enabled. However, depending on the specific use case and the hardware available,
+they may in fact slow down the execution of your script, because in both cases, there is an overhead cost involved:
+
+- **Parallel computation**: there is an overhead cost associated with initialising the separate threads.
+- **Numba JIT compilation**: there is an overhead cost associated with compiling the code to machine code. This happens only
+  the first time the function is called, but if the function is not called many times, this overhead cost can outweigh the
+  time saved on function execution.
+
+To experiment with what settings work best for your script, you can toggle these options on and off as follows:
+
+- **Disable parallel computation**: set the `parallel` user option to `False`:
+
+    .. code-block::
+
+         from rayflare.options import default_options
+         options = default_options
+         options.parallel = False
+
+
+  You can also control the number of workers which are spawned by setting the `n_jobs` user option. By default, it is
+  set to -1, which means all available CPUs are used, but it can also be set to a specific number.
+
+- **Disable Numba JIT compilation**: add the following at the start of your script - note that this must be before
+  any RayFlare classes/functions are imported:
+
+    .. code-block::
+
+         from numba import config
+         config.DISABLE_JIT = True

docs/Options/user_options.rst

@@ -2,7 +2,12 @@ Options
 =======
 .. _options:
 
-The function default_options is provided so not every
+.. toctree::
+   :maxdepth: 1
+
+   performance
+
+The function default_options (:literal:`from rayflare.options import default_options`) is provided so not every
 value has to be set manually by the user each time. All the user options, and their meanings, are listed below. As a general
 rule, values with units of distance (wavelength, depth spacing) are in m while angles are in radians.
 

docs/news.rst

@@ -11,12 +11,13 @@ To update to the latest version of RayFlare, run the following command in your t
 
    pip install rayflare --upgrade
 
-Version 2.0.0 released (2024-09-23)
+Version 2.0 released (2024-09-23)
 ------------------------------------
 **Highlights:**
 
 - **Speed improvements for all ray-tracing calculations**, by using numba jit (just-in-time) compilation of heavily-used
-  functions, and more efficient lookup of TMM values for integrated ray-tracing/TMM calculations.
+  functions, and more efficient lookup of TMM values for integrated ray-tracing/TMM calculations. The use of just-in-time
+  compilation can be turned on or off, see here: :ref:`Performance <performance>`.
 - **New analytical method for ray-tracing**, which is much faster than full ray-tracing and can be used with full
   accuracy for surfaces where the number of ray interactions is known exactly in advance (e.g. upright pyramids with
   opening angles between 45 and 54 degrees). Note that the calculation of absorption profiles in interface layers

examples/pvk_Si_analytical_RT.py

@@ -1,6 +1,15 @@
 from solcore.light_source import LightSource
 import numpy as np
-from rayflare.textures import regular_pyramids, planar_surface
+
+# the ray-tracer (from v2.0.0) will use numba's just-in-time (jit) compilation to make
+# heavily-used functions faster. However, the initial compilation also takes time, and thus if
+# these functions are not called many times, it is faster to disable the jit compilation
+# completely. The two lines below do this; this has to be done before importing any rayflare
+# functions or classes.
+from numba import config
+config.DISABLE_JIT = True
+
+from rayflare.textures import regular_pyramids
 import matplotlib.pyplot as plt
 
 from solcore import material
@@ -19,13 +28,12 @@
 GaAs = material("GaAs")()
 
 Pvk = material("Perovskite_CsBr_1p6eV")()
-glass = material("BK7")()
 
 n_rays = 2000
 
 d = 100e-6
 
-wavelengths = np.linspace(300, 1200, 50) * 1e-9
+wavelengths = np.linspace(300, 1000, 50) * 1e-9
 
 AM15G = LightSource(source_type='standard', version='AM1.5g', x=wavelengths,
                     output_units='photon_flux_per_m')
@@ -37,29 +45,30 @@
 options.nx = 20
 options.ny = 20
 options.theta_in = 0.1
-options.parallel = True
+options.parallel = False # for this example, with this choice of wavelengths, initializing the parallel
+# threads takes more time than it saves on executing the ray-tracing. For a more transparent structure
+# (at longer wavelengths in this case), this would no longer be true.
 options.randomize_surface = True
 options.I_thresh = 1e-3
 options.depth_spacing_bulk = 1e-8
 options.n_rays = n_rays
 options.maximum_passes = 100
 options.n_rays = n_rays
 options.pol = 'u'
+options.use_numba = False
 
-front_text = planar_surface(analytical=True)
-
-front_text_2 = regular_pyramids(52, True, 1,
+front_text = regular_pyramids(52, True, 1,
                                 interface_layers=[Layer(100e-9, MgF2), Layer(1000e-9, Pvk)],
                                 analytical=True)
 
 # note: not setting analytical = True for the last surface; rays will no longer be travelling in a
 # single direction after interacting with front_text_2, so this surface will not be treated analytically
-# anyway.
+# anyway, even if we set analytical = True.
 
 rear_text = regular_pyramids(52, False, 1)
 
-rt_str = rt_structure(textures=[front_text, front_text_2, rear_text], materials=[glass, Si],
-                      widths=[1000e-6, d], incidence=Air, transmission=Ag,
+rt_str = rt_structure(textures=[front_text, rear_text], materials=[Si],
+                      widths=[d], incidence=Air, transmission=Ag,
                       options=options, use_TMM=True, save_location='current',
                       overwrite=True)
 
@@ -69,17 +78,15 @@
 
 A_layer = result['A_per_layer']
 A_per_interface = result['A_per_interface']
-total = result['R'] + result['T'] + np.sum(A_layer, axis=1) + A_per_interface[1][:, 1]
-
+total = result['R'] + result['T'] + np.sum(A_layer, axis=1) + A_per_interface[0][:, 1]
 
 plt.figure()
 
 plt.plot(wavelengths * 1e9, result['R'], 'k-', label='R')
 plt.plot(wavelengths * 1e9, result['R0'], 'k-', alpha=0.5, label='R0')
 plt.plot(wavelengths * 1e9, result['T'], 'r-', label='T')
-plt.plot(wavelengths*1e9, A_layer[:,1], 'g-', label='Si')
-plt.plot(wavelengths*1e9, A_layer[:,0], 'b-', label='glass')
-plt.plot(wavelengths*1e9, A_per_interface[1][:,1], 'y-', label='Pvk')
+plt.plot(wavelengths*1e9, A_layer[:,0], 'b-', label='Si')
+plt.plot(wavelengths*1e9, A_per_interface[0][:,1], 'y-', label='Pvk')
 plt.plot(wavelengths * 1e9, total, 'k--', label='total', alpha=0.5)
 plt.xlabel("Wavelength (nm)")
 plt.ylabel("R/A/T")
@@ -92,6 +99,7 @@
 plt.plot(wavelengths * 1e9, np.max(result.ray_data.n_passes, axis=1), 'r-', label='Max. number of passes')
 plt.xlabel("Wavelength (nm)")
 plt.ylabel("Number of passes through a bulk material")
+plt.legend()
 plt.show()
 
 frac_reaching_max = [np.sum(x == options.maximum_passes) for x in result.ray_data.n_passes]

rayflare/ray_tracing/rt_common.py

@@ -10,10 +10,11 @@
 from math import atan2
 from random import random
 from copy import deepcopy
-from numba import jit
 from scipy.spatial import Delaunay
 from rayflare.utilities import process_pol
 
+from numba import jit
+
 unit_cell_N = np.array(
     [[0, -1, 0], [-1, 0, 0], [0, 1, 0], [1, 0, 0]]
 )  # surface normals: top, right, bottom, left
@@ -377,7 +378,7 @@ def decide_RT_TMM(ray, n0, n1, theta, N, side, rnd, lookuptable, d_theta):
 
     return side, A
 
-#@jit(nopython=True)
+@jit(nopython=True)
 def get_RT_data(theta, d_theta, R_data, T_data, Alayer_data, pol):
     # theta HAS to be positive here!
     angle_ind = int(np.floor(theta/d_theta)) # floor to avoid issues when theta = np.pi/2