3.3. WPM 3D with storing intermediate data

3D drawings are very nice, but sometimes we have not enough memory to store data. In this case we can use the following functions at Scalar_fields_XY scheme: WPM and BPM.

It is quite easy to execute them at vacuum, as only a XY initial field is required. However it is also possible to use a 3d refraction index structure. In order to not storing it, a function with the XYZ definition is required. However, the data are only extracted for each z position, and then no problems with storage are produced.

In my compute I have got the following results: - 128 x 128 x 128 = 600 ms - 256 x 256 x 256 = 5 seconds - 512 x 512 x 512 = 44 seconds (I could not resolve it in XYZ). - 1024 x 1024 x 1024 = 7 minutes (I could not resolve it in XYZ).

3.3.1. Modules

[1]:
from diffractio import np, plt, sp, um, mm, degrees
from diffractio.scalar_fields_XYZ import Scalar_field_XYZ

from diffractio.scalar_masks_XY import Scalar_mask_XY
from diffractio.scalar_sources_XY import Scalar_source_XY
from diffractio.scalar_masks_XYZ import Scalar_mask_XYZ

from diffractio.utils_math import get_k, ndgrid, nearest
from diffractio.scalar_fields_XY import WPM_schmidt_kernel

import time
import sys
number of processors: 8
total memory        : 7.5 Gb
available memory    : 50 %
max frequency       : 4000 GHz

3.4. Structure

[16]:
num_z=256

x=np.linspace(-50*um,50*um, 256)
y=np.linspace(-50*um,50*um, 256)
z=np.linspace(0,45*um,num_z)
z_obs=z[-1]
wavelength=1*um
[17]:
t0=Scalar_mask_XY(x,y,wavelength)
t0.circle(r0=(0,0), radius=45*um, angle=0)

u0=Scalar_source_XY(x,y,wavelength)
u0.gauss_beam(r0=(0,0), w0=45*um)

u_iter=u0*t0
[18]:
def fn(x,y,z,wavelegth):
    u=Scalar_mask_XYZ(x,y, z, wavelength)
    u.sphere(r0=(0,0,25), radius=20*um, refraction_index=2, angles=(0,0,0));
    index = u.n
    del u
    return np.squeeze(index)
[19]:
# Example of structure at a given k-position in z distance

refraction_index=fn(x,y,np.array([44,]),wavelength)
plt.imshow(refraction_index.real);
../../../_images/source_tutorial_scalar_XYZ_wpm_without_storing_7_0.png
[20]:
%%time
u_iter.WPM(fn, z_ini=z[0], z_end= z[-1], dz=(z[-1]-z[0])/num_z,
                            has_edges=True, matrix=False, verbose=True)

Time = 3.81 s, time/loop = 14.87 ms
CPU times: user 3.83 s, sys: 23.5 ms, total: 3.86 s
Wall time: 3.82 s
[21]:
u_iter.draw(has_colorbar='vertical', logarithm=True);
../../../_images/source_tutorial_scalar_XYZ_wpm_without_storing_9_0.png

We perform the same experiment, but storing data. Finally, a comparison between both procedures is performed.

[22]:
u_3d=Scalar_mask_XYZ(x,y, z, wavelength)
u_3d.n=fn(x,y,z,wavelength)
[23]:
u_3d.incident_field(u0=u0*t0)
[24]:
%%time
u_3d.clear_field()
u_3d.WPM(verbose=True, has_edges=True)
u_3d.draw_XZ(y0=0, logarithm=True)
Time = 6.58 s, time/loop = 25.71 ms
CPU times: user 6.73 s, sys: 58.2 ms, total: 6.79 s
Wall time: 6.71 s
[24]:
<matplotlib.image.AxesImage at 0x7f1416963160>
../../../_images/source_tutorial_scalar_XYZ_wpm_without_storing_13_2.png
[25]:
u_3d.draw_XY(z0=50, logarithm=True, normalize='',has_colorbar='vertical')
../../../_images/source_tutorial_scalar_XYZ_wpm_without_storing_14_0.png

3.4.1. Differences between both algorithms

[26]:
I_3d_final=np.abs(u_3d.u[:,:,-1])**2
I_iter=np.abs(u_iter.u)**2

u_diffs=Scalar_mask_XY(x,y,wavelength)
u_diffs.u= (I_iter-I_3d_final)/I_iter.max()

u_diffs.draw(has_colorbar='vertical', logarithm=False)
print("Differences = {:2.2f} %".format((100*np.abs((I_iter-I_3d_final)/I_iter.max()).mean())))
Differences = 0.01 %
../../../_images/source_tutorial_scalar_XYZ_wpm_without_storing_16_1.png
[ ]: