From 9661fb0cfd543fbd545fa3863250bb1503d267c5 Mon Sep 17 00:00:00 2001 From: Logan Bishop-Van Horn Date: Wed, 4 Oct 2023 11:00:04 -0700 Subject: [PATCH] More memory efficient biot_savart_2d (#107) * More memory efficient biot_savart_2d --- docs/notebooks/field-sources.ipynb | 10 +- docs/notebooks/logo.ipynb | 6 +- docs/notebooks/polygons.ipynb | 20 ++-- docs/notebooks/quickstart.ipynb | 10 +- docs/notebooks/scanning-squid.ipynb | 12 +-- docs/notebooks/terminal-currents.ipynb | 6 +- superscreen/distance.py | 46 +++------ superscreen/sources/current.py | 134 ++++++++++++++++++++----- superscreen/test/test_distance.py | 9 -- superscreen/version.py | 2 +- 10 files changed, 157 insertions(+), 98 deletions(-) diff --git a/docs/notebooks/field-sources.ipynb b/docs/notebooks/field-sources.ipynb index 90ee49a..71ccfee 100644 --- a/docs/notebooks/field-sources.ipynb +++ b/docs/notebooks/field-sources.ipynb @@ -87,8 +87,8 @@ { "data": { "text/html": [ - "
SoftwareVersion
SuperScreen0.10.2
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", - "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Tue Sep 19 16:13:07 2023 PDT
" + "
SoftwareVersion
SuperScreen0.10.3
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", + "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Wed Oct 04 10:31:44 2023 PDT
" ], "text/plain": [ "" @@ -142,9 +142,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "x = 0.7650684631746568\n", - "y = 0.7478751704881408\n", - "z = 0.18189859902136085\n", + "x = 0.3396785462775239\n", + "y = 0.5560346157741703\n", + "z = 0.47978056547366055\n", "field(x, y, z) = 5.0\n" ] } diff --git a/docs/notebooks/logo.ipynb b/docs/notebooks/logo.ipynb index ce5c14b..6014cd0 100644 --- a/docs/notebooks/logo.ipynb +++ b/docs/notebooks/logo.ipynb @@ -47,7 +47,7 @@ "\n", "import superscreen as sc\n", "\n", - "SAVE = False" + "SAVE = True" ] }, { @@ -59,8 +59,8 @@ { "data": { "text/html": [ - "
SoftwareVersion
SuperScreen0.10.2
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", - "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Tue Sep 19 16:13:38 2023 PDT
" + "
SoftwareVersion
SuperScreen0.10.3
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", + "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Wed Oct 04 10:32:27 2023 PDT
" ], "text/plain": [ "" diff --git a/docs/notebooks/polygons.ipynb b/docs/notebooks/polygons.ipynb index 88fb58a..7d24ad3 100644 --- a/docs/notebooks/polygons.ipynb +++ b/docs/notebooks/polygons.ipynb @@ -60,8 +60,8 @@ { "data": { "text/html": [ - "
SoftwareVersion
SuperScreen0.10.2
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", - "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Tue Sep 19 16:14:10 2023 PDT
" + "
SoftwareVersion
SuperScreen0.10.3
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", + "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Wed Oct 04 10:31:31 2023 PDT
" ], "text/plain": [ "" @@ -117,7 +117,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 5, @@ -209,7 +209,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 7, @@ -235,7 +235,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 8, @@ -264,7 +264,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 9, @@ -290,7 +290,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 10, @@ -316,7 +316,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 11, @@ -378,7 +378,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 13, @@ -402,7 +402,7 @@ "" ], "text/plain": [ - "" + "" ] }, "execution_count": 14, diff --git a/docs/notebooks/quickstart.ipynb b/docs/notebooks/quickstart.ipynb index 5358cc7..b074c50 100644 --- a/docs/notebooks/quickstart.ipynb +++ b/docs/notebooks/quickstart.ipynb @@ -55,8 +55,8 @@ { "data": { "text/html": [ - "
SoftwareVersion
SuperScreen0.10.2
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", - "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Tue Sep 19 16:13:03 2023 PDT
" + "
SoftwareVersion
SuperScreen0.10.3
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", + "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Wed Oct 04 10:28:14 2023 PDT
" ], "text/plain": [ "" @@ -1105,7 +1105,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Holes: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:01<00:00, 1.92it/s]" + "Holes: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 2.55it/s]" ] }, { @@ -1193,7 +1193,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Holes: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 2.14it/s]\n" + "Holes: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 2.77it/s]\n" ] } ], @@ -1430,7 +1430,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Solver iterations: 100%|█████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [00:02<00:00, 3.72it/s]\n" + "Solver iterations: 100%|█████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [00:02<00:00, 4.06it/s]\n" ] } ], diff --git a/docs/notebooks/scanning-squid.ipynb b/docs/notebooks/scanning-squid.ipynb index 4b43acc..d9f4fb3 100644 --- a/docs/notebooks/scanning-squid.ipynb +++ b/docs/notebooks/scanning-squid.ipynb @@ -97,8 +97,8 @@ { "data": { "text/html": [ - "
SoftwareVersion
SuperScreen0.10.2
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", - "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Tue Sep 19 16:14:27 2023 PDT
" + "
SoftwareVersion
SuperScreen0.10.3
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", + "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Wed Oct 04 10:29:53 2023 PDT
" ], "text/plain": [ "" @@ -160,10 +160,10 @@ "name": "stderr", "output_type": "stream", "text": [ - "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:01<00:00, 2.77it/s]\n", - "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:01<00:00, 2.82it/s]\n", - "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00, 1.36it/s]\n", - "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00, 1.52it/s]\n" + "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:01<00:00, 2.65it/s]\n", + "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:01<00:00, 2.81it/s]\n", + "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00, 1.41it/s]\n", + "Solver iterations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:03<00:00, 1.61it/s]\n" ] } ], diff --git a/docs/notebooks/terminal-currents.ipynb b/docs/notebooks/terminal-currents.ipynb index 635b0b3..84d894f 100644 --- a/docs/notebooks/terminal-currents.ipynb +++ b/docs/notebooks/terminal-currents.ipynb @@ -157,8 +157,8 @@ { "data": { "text/html": [ - "
SoftwareVersion
SuperScreen0.10.2
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", - "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Tue Sep 19 16:14:52 2023 PDT
" + "
SoftwareVersion
SuperScreen0.10.3
Numpy1.23.3
Numba0.57.0
SciPy1.9.1
matplotlib3.6.0
IPython8.5.0
Python3.9.13 | packaged by conda-forge | (main, May 27 2022, 17:01:00) \n", + "[Clang 13.0.1 ]
OSposix [darwin]
Number of CPUsPhysical: 10, Logical: 10
BLAS InfoOPENBLAS
Wed Oct 04 10:29:14 2023 PDT
" ], "text/plain": [ "" @@ -480,7 +480,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Holes: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 2.91it/s]\n" + "Holes: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 3.01it/s]\n" ] } ], diff --git a/superscreen/distance.py b/superscreen/distance.py index a9a1109..90bb1c5 100644 --- a/superscreen/distance.py +++ b/superscreen/distance.py @@ -2,33 +2,15 @@ import numpy as np -@numba.njit(fastmath=True, parallel=True) -def pairwise_difference(xA: np.ndarray, xB: np.ndarray): - """Pairwise different between two 1D arrays. - - This is equivalent to ``numpy.subtract.outer(XA, XB)``. - - Args: - xA: A shape (n,) array - xB: A shep (m,) array - - Returns: - A shape (n, m) array of pairwise differences. - """ - out = np.empty((xA.shape[0], xB.shape[0]), dtype=xA.dtype) - for i in numba.prange(xA.shape[0]): - for j in range(xB.shape[0]): - out[i, j] = xA[i] - xB[j] - return out - - @numba.njit(fastmath=True, parallel=True) def sqeuclidean_distance_2d(XA: np.ndarray, XB: np.ndarray): """Squared Euclidean pointwise distance between two 2D arrays.""" out = np.empty((XA.shape[0], XB.shape[0]), dtype=XA.dtype) for i in numba.prange(XA.shape[0]): for j in range(XB.shape[0]): - out[i, j] = (XA[i, 0] - XB[j, 0]) ** 2 + (XA[i, 1] - XB[j, 1]) ** 2 + dx = XA[i, 0] - XB[j, 0] + dy = XA[i, 1] - XB[j, 1] + out[i, j] = dx * dx + dy * dy return out @@ -38,11 +20,10 @@ def sqeuclidean_distance_3d(XA: np.ndarray, XB: np.ndarray): out = np.empty((XA.shape[0], XB.shape[0]), dtype=XA.dtype) for i in numba.prange(XA.shape[0]): for j in range(XB.shape[0]): - out[i, j] = ( - (XA[i, 0] - XB[j, 0]) ** 2 - + (XA[i, 1] - XB[j, 1]) ** 2 - + (XA[i, 2] - XB[j, 2]) ** 2 - ) + dx = XA[i, 0] - XB[j, 0] + dy = XA[i, 1] - XB[j, 1] + dz = XA[i, 2] - XB[j, 2] + out[i, j] = dx * dx + dy * dy + dz * dz return out @@ -52,7 +33,9 @@ def euclidean_distance_2d(XA: np.ndarray, XB: np.ndarray): out = np.empty((XA.shape[0], XB.shape[0]), dtype=XA.dtype) for i in numba.prange(XA.shape[0]): for j in range(XB.shape[0]): - out[i, j] = np.sqrt((XA[i, 0] - XB[j, 0]) ** 2 + (XA[i, 1] - XB[j, 1]) ** 2) + dx = XA[i, 0] - XB[j, 0] + dy = XA[i, 1] - XB[j, 1] + out[i, j] = np.sqrt(dx * dx + dy * dy) return out @@ -62,11 +45,10 @@ def euclidean_distance_3d(XA: np.ndarray, XB: np.ndarray): out = np.empty((XA.shape[0], XB.shape[0]), dtype=XA.dtype) for i in numba.prange(XA.shape[0]): for j in range(XB.shape[0]): - out[i, j] = np.sqrt( - (XA[i, 0] - XB[j, 0]) ** 2 - + (XA[i, 1] - XB[j, 1]) ** 2 - + (XA[i, 2] - XB[j, 2]) ** 2 - ) + dx = XA[i, 0] - XB[j, 0] + dy = XA[i, 1] - XB[j, 1] + dz = XA[i, 2] - XB[j, 2] + out[i, j] = np.sqrt(dx * dx + dy * dy + dz * dz) return out diff --git a/superscreen/sources/current.py b/superscreen/sources/current.py index eaccc21..27dac23 100644 --- a/superscreen/sources/current.py +++ b/superscreen/sources/current.py @@ -1,15 +1,115 @@ from typing import Optional, Union +import numba import numpy as np from scipy.constants import mu_0 from scipy.spatial import Delaunay from ..device.utils import vertex_areas -from ..distance import pairwise_difference from ..parameter import Parameter from ..units import ureg +@numba.njit(fastmath=True, parallel=True) +def _biot_savart_2d_z( + eval_positions: np.ndarray, + positions: np.ndarray, + current_densities: np.ndarray, + areas: np.ndarray, +) -> np.ndarray: + """Returns the z-component of the magnetic field (in tesla) + from a given sheet current density distribution. + + Args: + eval_positions: The ``(x, y, z)`` coordinates at which to evaluate the field in meters + positions: The ``(x, y, z)`` coordinates of the sheet current in meters + current_densities: The sheet current density in amps / meter + areas: The effective areas of the sheet current mesh in meters**2 + + Returns: + The z-component of the magnetic field in tesla + """ + assert eval_positions.ndim == 2 + assert eval_positions.shape[1] == 3 + assert positions.ndim == 2 + assert positions.shape[0] == areas.shape[0] == current_densities.shape[0] + assert positions.shape[1] == 3 + + Jx = current_densities[:, 0] + Jy = current_densities[:, 1] + Bz_out = np.empty(len(eval_positions), dtype=float) + + for i in numba.prange(eval_positions.shape[0]): + Jx_dy = 0.0 + Jy_dx = 0.0 + for k in range(positions.shape[0]): + dx = eval_positions[i, 0] - positions[k, 0] + dy = eval_positions[i, 1] - positions[k, 1] + dz = eval_positions[i, 2] - positions[k, 2] + pref = ( + (mu_0 / (4 * np.pi)) + * areas[k] + * (dx * dx + dy * dy + dz * dz) ** (-3 / 2) + ) + Jx_dy += pref * Jx[k] * dy + Jy_dx += pref * Jy[k] * dx + Bz_out[i] = Jx_dy - Jy_dx + return Bz_out + + +@numba.njit(fastmath=True, parallel=True) +def _biot_savart_2d_vector( + eval_positions: np.ndarray, + positions: np.ndarray, + current_densities: np.ndarray, + areas: np.ndarray, +) -> np.ndarray: + """Returns the vector magnetic field (in tesla) + from a given sheet current density distribution. + + Args: + eval_positions: The ``(x, y, z)`` coordinates at which to evaluate the field in meters + positions: The ``(x, y, z)`` coordinates of the sheet current in meters + current_densities: The sheet current density in amps / meter + areas: The effective areas of the sheet current mesh in meters**2 + + Returns: + The vector magnetic field in tesla + """ + assert eval_positions.ndim == 2 + assert eval_positions.shape[1] == 3 + assert positions.ndim == 2 + assert positions.shape[0] == areas.shape[0] == current_densities.shape[0] + assert positions.shape[1] == 3 + + Jx = current_densities[:, 0] + Jy = current_densities[:, 1] + B_out = np.empty((len(eval_positions), 3), dtype=float) + + for i in numba.prange(eval_positions.shape[0]): + Jx_dy = 0.0 + Jy_dx = 0.0 + Jx_dz = 0.0 + Jy_dz = 0.0 + for k in range(positions.shape[0]): + dx = eval_positions[i, 0] - positions[k, 0] + dy = eval_positions[i, 1] - positions[k, 1] + dz = eval_positions[i, 2] - positions[k, 2] + pref = ( + (mu_0 / (4 * np.pi)) + * areas[k] + * (dx * dx + dy * dy + dz * dz) ** (-3 / 2) + ) + Jx_dy += pref * Jx[k] * dy + Jy_dx += pref * Jy[k] * dx + Jx_dz += pref * Jx[k] * dz + Jy_dz += pref * Jy[k] * dz + B_out[i, 0] = Jy_dz + B_out[i, 1] = -Jx_dz + B_out[i, 2] = Jx_dy - Jy_dx + return B_out + + def biot_savart_2d( x: Union[float, np.ndarray], y: Union[float, np.ndarray], @@ -76,38 +176,24 @@ def biot_savart_2d( x, y, z = np.atleast_1d(x, y, z) if z.shape[0] == 1: z = z * np.ones_like(x) - x = x * to_meter - y = y * to_meter - z = z * to_meter + eval_positions = np.array([x, y, z]).T * to_meter positions, current_densities = np.atleast_2d(positions, current_densities) - positions = positions * to_meter - z0 = z0 * to_meter current_densities = current_densities * to_amp_per_meter - # Calculate the pairwise distance between the current sheet and evaluation - # points for each axis. - x0, y0 = positions[:, 0], positions[:, 1] - Jx, Jy = current_densities[:, 0], current_densities[:, 1] - dx = pairwise_difference(x, x0) - dy = pairwise_difference(y, y0) - dz = pairwise_difference(z, z0 * np.ones_like(x0)) + positions = positions * to_meter + z0 = z0 * np.ones(len(positions)) * to_meter if areas is None: # Triangulate the current sheet to assign an effective area to each vertex. triangles = Delaunay(positions).simplices areas = vertex_areas(positions, triangles) else: areas = areas * to_meter**2 + positions = np.concatenate([positions, z0[:, np.newaxis]], axis=1) # Evaluate the Biot-Savart integral. - pref = (mu_0 / (4 * np.pi)) * areas * (dx**2 + dy**2 + dz**2) ** (-3 / 2) - Jx_dy = np.einsum("ij, ij, j -> i", pref, dy, Jx) - Jy_dx = np.einsum("ij, ij, j -> i", pref, dx, Jy) - Bz = Jx_dy - Jy_dx - if not vector: - return Bz - Jy_dz = np.einsum("ij, ij, j -> i", pref, dz, Jy) - Jx_dz = np.einsum("ij, ij, j -> i", pref, dz, Jx) - Bx = Jy_dz - By = -Jx_dz - return np.stack([Bx, By, Bz], axis=1) + if vector: + B = _biot_savart_2d_vector(eval_positions, positions, current_densities, areas) + else: + B = _biot_savart_2d_z(eval_positions, positions, current_densities, areas) + return B def SheetCurrentField( diff --git a/superscreen/test/test_distance.py b/superscreen/test/test_distance.py index a3493f0..ffd1d3b 100644 --- a/superscreen/test/test_distance.py +++ b/superscreen/test/test_distance.py @@ -35,12 +35,3 @@ def test_cdist(metric, dtype, ndim): dist_sc = sc.distance.cdist(XA, XB, metric=metric) dist_sp = scipy.spatial.distance.cdist(XA, XB, metric=metric) assert np.allclose(dist_sc, dist_sp) - - -@pytest.mark.parametrize("dtype", ("float64", "float32")) -def test_pairwise_difference(dtype): - XA = np.random.random(400).astype(dtype) - XB = np.random.random(777).astype(dtype) - diff_sc = sc.distance.pairwise_difference(XA, XB) - diff_np = np.subtract.outer(XA, XB) - assert np.allclose(diff_sc, diff_np) diff --git a/superscreen/version.py b/superscreen/version.py index ae69f98..11ef801 100644 --- a/superscreen/version.py +++ b/superscreen/version.py @@ -1,2 +1,2 @@ -__version_info__ = (0, 10, 2) +__version_info__ = (0, 10, 3) __version__ = ".".join(map(str, __version_info__))