From 1dd34dd2663f2a52ccc1848967be178d38356383 Mon Sep 17 00:00:00 2001 From: Nick Karpowicz <101142784+NickKarpowicz@users.noreply.github.com> Date: Thu, 12 Feb 2026 13:01:33 +0100 Subject: [PATCH 1/5] use std::aligned_alloc when allocating CPU device memory --- Source/Devices/LWEActiveDeviceCPU.h | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Source/Devices/LWEActiveDeviceCPU.h b/Source/Devices/LWEActiveDeviceCPU.h index 72472a48..60362f88 100644 --- a/Source/Devices/LWEActiveDeviceCPU.h +++ b/Source/Devices/LWEActiveDeviceCPU.h @@ -246,18 +246,23 @@ class CPUDevice : public LWEDevice { } #endif int deviceCalloc(void** ptr, const size_t N, const size_t elementSize) override { - *ptr = calloc(N, elementSize); - return static_cast(*ptr == nullptr); + const size_t request = N*elementSize; + const size_t standard_alignment = 2*sizeof(deviceFP); + const size_t bytes = (request % standard_alignment) ? standard_alignment * (request / standard_alignment) + standard_alignment : request; + *ptr = std::aligned_alloc(standard_alignment, bytes); + if(*ptr == nullptr) return 1; + std::memset(*ptr, 0, bytes); + return 0; } void deviceMemset(void* ptr, int value, size_t count) override { - memset(ptr, value, count); + std::memset(ptr, value, count); } void deviceMemcpyImplementation(void* dst, const void* src, size_t count, copyType kind) override { std::memcpy(dst, src, count); } void deviceFree(void* block) override { - free(block); + std::free(block); } inline bool isTheCanaryPixelNaN(const deviceFP* canaryPointer) { From 0c50de12fff2c8c8c27e06f04e094dd27a57386f Mon Sep 17 00:00:00 2001 From: Nick Karpowicz <101142784+NickKarpowicz@users.noreply.github.com> Date: Mon, 16 Feb 2026 10:05:45 +0100 Subject: [PATCH 2/5] fix addPulse, fixes #95 --- .../LightwaveExplorerDocumentation.py | 132 +++++++++++++----- Source/LightwaveExplorerCore.cu | 17 +-- 2 files changed, 106 insertions(+), 43 deletions(-) diff --git a/Documentation/LightwaveExplorerDocumentation.py b/Documentation/LightwaveExplorerDocumentation.py index 6909dd98..ccffaae7 100644 --- a/Documentation/LightwaveExplorerDocumentation.py +++ b/Documentation/LightwaveExplorerDocumentation.py @@ -111,7 +111,7 @@ def _(mo): @app.cell(hide_code=True) def _(mo): - #Slider to control the UV resonance + # Slider to control the UV resonance UV = mo.ui.slider(start=0.5, stop=2.0, step=0.1, value=1, label="UV factor") UV return (UV,) @@ -119,7 +119,7 @@ def _(mo): @app.cell(hide_code=True) def _(mo): - #Slider to control the IR resonance + # Slider to control the IR resonance IR = mo.ui.slider(start=0.5, stop=2.0, step=0.1, value=1, label="IR factor") IR return (IR,) @@ -127,45 +127,104 @@ def _(mo): @app.cell(hide_code=True) def _(IR, UV, mo): - #importing some modules I'll want - import LightwaveExplorer as lwe + # importing some modules I'll want import io - import numpy as np + + import LightwaveExplorer as lwe import matplotlib.pyplot as plt + import numpy as np from matplotlib import rcParams - rcParams['font.family'] = 'sans-serif' - rcParams['font.sans-serif'] = ['Helvetica', 'Arial', 'Verdana', 'DejaVu Sans', 'Liberation Sans', 'Bitstream Vera Sans', 'sans-serif'] + rcParams["font.family"] = "sans-serif" + rcParams["font.sans-serif"] = [ + "Helvetica", + "Arial", + "Verdana", + "DejaVu Sans", + "Liberation Sans", + "Bitstream Vera Sans", + "sans-serif", + ] def showmo(): """ Helper function to plot as an svg and have it display in marimo in vector form """ svg_buffer = io.StringIO() - plt.savefig(svg_buffer, format='svg') + plt.savefig(svg_buffer, format="svg") svg_buffer.seek(0) svg_data = svg_buffer.getvalue() return mo.Html(svg_data) - #first we'll make a wavelength grid to work with - l = np.linspace(0.3,3,1024) - - #next we'll need Sellmeier coefficients, these are for barium fluoride, H. H. Li., J. Phys. Chem. Ref. Data 9, 161-289 (1980) - a = np.array([1.33973,0,0.81070,-0.010130,0,0.19652,-892.22,0,4.52469,-2896.6,0,0,1,0,0,0,0,0,0,0,0,0]) - - #we can get the refractive index for the wavelengths we put in the grid by calling the sellmeier() function - #from the lightwaveExplorer module, with the equationType set to 0. + # first we'll make a wavelength grid to work with + l = np.linspace(0.3, 3, 1024) + + # next we'll need Sellmeier coefficients, these are for barium fluoride, H. H. Li., J. Phys. Chem. Ref. Data 9, 161-289 (1980) + a = np.array( + [ + 1.33973, + 0, + 0.81070, + -0.010130, + 0, + 0.19652, + -892.22, + 0, + 4.52469, + -2896.6, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + ] + ) + + # we can get the refractive index for the wavelengths we put in the grid by calling the sellmeier() function + # from the lightwaveExplorer module, with the equationType set to 0. n = lwe.sellmeier(l, a, 0) - #let's make it so we can adjust the resonances of the oscillators and see how it affects the index - - fig,ax = plt.subplots(1,1, figsize=(5,4)) - a2 = np.array([1.33973,0,0.81070,-0.010130,0,0.19652,-892.22,0,4.52469,-2896.6,0,0,1,0,0,0,0,0,0,0,0,0]) + # let's make it so we can adjust the resonances of the oscillators and see how it affects the index + + fig, ax = plt.subplots(1, 1, figsize=(5, 4)) + a2 = np.array( + [ + 1.33973, + 0, + 0.81070, + -0.010130, + 0, + 0.19652, + -892.22, + 0, + 4.52469, + -2896.6, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + ] + ) a2[3] *= UV.value a2[6] *= IR.value n2 = lwe.sellmeier(l, a2, 0) - ax.plot(l,np.real(n),label="original",color="blue") - ax.plot(l,np.real(n2), label = "modified", color = "magenta") + ax.plot(l, np.real(n), label="original", color="blue") + ax.plot(l, np.real(n2), label="modified", color="magenta") ax.set_xlabel("Wavelength (microns)") ax.set_ylabel("n") ax.legend() @@ -549,16 +608,23 @@ def _(mo): def _(Ns, bandwidth, cep, f0, lwe, np, phi2, phi3, plt, showmo, tau): dt = 0.25e-15 Tgrid = 80e-15 - Nt = int(2*Tgrid/dt) - t = np.linspace(-Tgrid,Tgrid,Nt) - w = 2*np.pi*np.fft.fftfreq(Nt,dt) - ws = (w-2*np.pi*f0.value*1e12) - phi = np.pi*cep.value + (-Tgrid + tau.value*1e-15)*w + 0.5*phi2.value*1e-30*ws**2 + (1.0/6)*phi3.value*1e-45*ws**3 - Ew = np.exp(-ws**Ns.value/(2*np.pi*bandwidth.value*1e12)**Ns.value -1.0j * phi) - Ew[w<0]=0 + Nt = int(2 * Tgrid / dt) + t = np.linspace(-Tgrid, Tgrid, Nt) + w = 2 * np.pi * np.fft.fftfreq(Nt, dt) + ws = w - 2 * np.pi * f0.value * 1e12 + phi = ( + np.pi * cep.value + + (-Tgrid + tau.value * 1e-15) * w + + 0.5 * phi2.value * 1e-30 * ws**2 + + (1.0 / 6) * phi3.value * 1e-45 * ws**3 + ) + Ew = np.exp( + -(ws**Ns.value) / (2 * np.pi * bandwidth.value * 1e12) ** Ns.value - 1.0j * phi + ) + Ew[w < 0] = 0 Et = lwe.norma(np.fft.ifft(Ew)) - fwhm_value = 1e15*lwe.fwhm(t,np.abs(Et)**2) - plt.plot(1e15*t,np.real(Et)) + fwhm_value = 1e15 * lwe.fwhm(t, np.abs(Et) ** 2) + plt.plot(1e15 * t, np.real(Et)) plt.xlabel("Time (fs)") plt.title(f"FWHM duration: {fwhm_value: .2f} fs") showmo() @@ -803,7 +869,7 @@ def _(mo): Note that compared to the interface, you have an extra beam angle and beam position. These are only relevant to 3D calculations and gives the angle and offset in the y-direction. (you can actually specify them on the interface, too, by writing both angles in the text box separated by a semicolon in the box for x0 and the NC angle, e.g. 10;20 will give an x-direction value of 10 and a y-direction value of 20). - A call to addPulse() will initialize the electric field grid if it hasn't been done so already. + A call to addPulse() will initialize the electric field grid if it hasn't been done so already, with the added pulse only. #### linear(materialIndex, crystalTheta, crystalPhi, crystalThickness, propagationStep) This function will propagate *linearly* through the medium determined by materialIndex. The parameter propagationStep is only necessary when in cylindrical symmetry mode. This is because the propagation is carried out analytically in 3D and 2D Cartesian modes since the linear propagator is diagonal in those bases. In cylindrical symmetry mode, it is handled as a numerical propagation with all nonlinearities disabled, which is why specifying a step size is necessary. This is also why this function is significantly faster in 3D and 2D Cartesian coordinates. @@ -1003,7 +1069,7 @@ def _(mo): @app.cell(hide_code=True) def _(mo): mo.md(r""" - + """) return diff --git a/Source/LightwaveExplorerCore.cu b/Source/LightwaveExplorerCore.cu index b43fa9a2..435d4522 100644 --- a/Source/LightwaveExplorerCore.cu +++ b/Source/LightwaveExplorerCore.cu @@ -2444,7 +2444,6 @@ namespace hostFunctions{ copyType::ToDevice); Pulse p; - p = sCPU->pulse1; p.energy = parameters[0]; p.frequency = 1e12 * parameters[1]; p.bandwidth = 1e12 * parameters[2]; @@ -2455,15 +2454,13 @@ namespace hostFunctions{ p.tod = 1e-45 * parameters[7]; p.phaseMaterial = (int)parameters[8]; p.phaseMaterialThickness = 1e-6 * parameters[9]; - p.beam_spec.relevant_modes = 1; - p.beam_spec.relevant_expansion = 1; - p.beam_spec.weight[0] = 1.0; - p.beam_spec.waist[0][0] = 1e-6 * parameters[10]; - p.beam_spec.x_offset[0][0] = 1e-6 * parameters[11]; - p.beam_spec.y_offset[0][0] = 1e-6 * parameters[12]; - p.beam_spec.z_offset[0][0] = 1e-6 * parameters[13]; - p.beam_spec.angle_x[0][0] = deg2Rad() * parameters[14]; - p.beam_spec.angle_y[0][0] = deg2Rad() * parameters[15]; + p.beam_spec = BeamSpecification{}; + p.beamwaist = 1e-6 * parameters[10]; + p.x_offset = 1e-6 * parameters[11]; + p.y_offset = 1e-6 * parameters[12]; + p.z_offset = 1e-6 *parameters[13]; + p.angle_x_offset = deg2Rad() * parameters[14]; + p.angle_y_offset = deg2Rad() * parameters[15]; p.polarizationAngle = deg2Rad() * parameters[16]; p.circularity = parameters[17]; (*sCPU).materialIndex = (int)parameters[18]; From 0fb903b6c17416e1f913dd077cb52daaed98cb89 Mon Sep 17 00:00:00 2001 From: Nick Karpowicz <101142784+NickKarpowicz@users.noreply.github.com> Date: Mon, 16 Feb 2026 10:41:21 +0100 Subject: [PATCH 3/5] fix return type of BeamSpecification operator= --- Source/DataStructures.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/DataStructures.hpp b/Source/DataStructures.hpp index 3beffe24..bd69bf82 100644 --- a/Source/DataStructures.hpp +++ b/Source/DataStructures.hpp @@ -768,7 +768,7 @@ struct BeamSpecification { } template - BeamSpecification& operator=(const BeamSpecification& other){ + void operator=(const BeamSpecification& other){ basis = other.basis; relevant_modes = other.relevant_modes; relevant_expansion = other.relevant_expansion; From a85a81f4de7b3e2c725edb438ebca2c632b97e4b Mon Sep 17 00:00:00 2001 From: Nick Karpowicz <101142784+NickKarpowicz@users.noreply.github.com> Date: Mon, 16 Feb 2026 11:43:10 +0100 Subject: [PATCH 4/5] include for std;:aligned_alloc --- Source/Devices/LWEActiveDeviceCPU.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Devices/LWEActiveDeviceCPU.h b/Source/Devices/LWEActiveDeviceCPU.h index 60362f88..88bedddf 100644 --- a/Source/Devices/LWEActiveDeviceCPU.h +++ b/Source/Devices/LWEActiveDeviceCPU.h @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #ifdef __APPLE__ From 54476516baa1dc813dfb0dcaf152731b7840b16a Mon Sep 17 00:00:00 2001 From: Nick Karpowicz <101142784+NickKarpowicz@users.noreply.github.com> Date: Mon, 16 Feb 2026 12:32:41 +0100 Subject: [PATCH 5/5] don't use std::aligned_alloc on windows because it doesn't support it --- Source/Devices/LWEActiveDeviceCPU.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Source/Devices/LWEActiveDeviceCPU.h b/Source/Devices/LWEActiveDeviceCPU.h index 88bedddf..7491448e 100644 --- a/Source/Devices/LWEActiveDeviceCPU.h +++ b/Source/Devices/LWEActiveDeviceCPU.h @@ -5,7 +5,6 @@ #include #include #include -#include #include #include #ifdef __APPLE__ @@ -250,7 +249,11 @@ class CPUDevice : public LWEDevice { const size_t request = N*elementSize; const size_t standard_alignment = 2*sizeof(deviceFP); const size_t bytes = (request % standard_alignment) ? standard_alignment * (request / standard_alignment) + standard_alignment : request; - *ptr = std::aligned_alloc(standard_alignment, bytes); + #ifndef _WIN32 + *ptr = std::aligned_alloc(standard_alignment, bytes); + #else + *ptr = std::malloc(bytes); + #endif if(*ptr == nullptr) return 1; std::memset(*ptr, 0, bytes); return 0;