Python CPU和GPU加速TDEs，每秒超过1亿个TDEs！

2022年12月的注释：这里的代码运行良好，我计划继续进行一些小错误修复，以保持当前的功能。但我就此不再对项目进行任何改进。我很乐意 聊天或电子邮件 帮助任何想深入研究的人。对于最大的潜在改进，请参阅此问题： https://github.com/tbenthompson/cutde/issues/23

cutde：启用CUDA、OpenCL和C++的全空间和半空间三角形滑移元素（TDEs），基准测试每秒130亿个TDEs。 cutde 是Nikhoo和Walter 2015的原始MATLAB代码的翻译和优化。除了基本的成对TDE操作（位移和应变）外，cutde 还包括

• 所有成对矩阵构建函数。
• 无矩阵函数，适用于低内存使用设置。
• 分块函数，特别适用于FMM或分层矩阵设置。
• 自适应交叉逼近（ACA）实现，用于构建分层矩阵。

下面是基本使用和安装说明。对于更实际的用法示例，请查看 这些BIE教程中的TDE示例。您将找到使用上述所有变体的示例。

• Python CPU和GPU加速TDEs，每秒超过1亿个TDEs！
• 使用文档
  • 简单的成对TDEs
  • 我想得到应力。
  • 所有成对交互矩阵
  • 无矩阵所有成对交互
  • 分块交互矩阵
  • 自适应交叉逼近（ACA）
• 安装
  • GPU安装
    • PyCUDA
    • Mac OS X
    • Ubuntu + PyOpenCL/PoCL
    • Ubuntu + PyOpenCL与系统驱动程序
    • Windows
    • 其他
    • 为什么我不能使用Apple CPU OpenCL？
• 开发
  • 架构

import matplotlib.pyplot as plt
import numpy as np

import cutde.fullspace as FS

xs = np.linspace(-2, 2, 200)
ys = np.linspace(-2, 2, 200)
obsx, obsy = np.meshgrid(xs, ys)
pts = np.array([obsx, obsy, 0 * obsy]).reshape((3, -1)).T.copy()

fault_pts = np.array([[-1, 0, 0], [1, 0, 0], [1, 0, -1], [-1, 0, -1]])
fault_tris = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64)

disp_mat = FS.disp_matrix(obs_pts=pts, tris=fault_pts[fault_tris], nu=0.25)

slip = np.array([[1, 0, 0], [1, 0, 0]])
disp = disp_mat.reshape((-1, 6)).dot(slip.flatten())

disp_grid = disp.reshape((*obsx.shape, 3))

plt.figure(figsize=(5, 5), dpi=300)
cntf = plt.contourf(obsx, obsy, disp_grid[:, :, 0], levels=21)
plt.contour(
    obsx,
    obsy,
    disp_grid[:, :, 0],
    colors="k",
    linestyles="-",
    linewidths=0.5,
    levels=21,
)
plt.colorbar(cntf)
plt.title("$u_x$")
plt.tight_layout()
plt.savefig("docs/example.png", bbox_inches="tight")

使用文档

简单的成对TDEs

为观测点/源元素对计算TDEs非常简单

from cutde.fullspace import disp, strain

disp = disp(obs_pts, src_tris, slips, nu)
strain = strain(obs_pts, src_tris, slips, nu)

将cutde.fullspace替换为cutde.halfspace以使用半空间TDEs。

参数：

• obs_pts 是一个形状为 (N, 3) 的 np.array。
• src_tris 是一个形状为 (N, 3, 3) 的 np.array，其中第二个维度对应于每个顶点，第三个维度对应于这些顶点的坐标。
• slips 是一个形状为 (N, 3) 的 np.array，其中 slips[:,0] 是剪切滑移分量，分量1是倾滑，分量2是拉伸/张开分量。
• 最后一个参数，nu，是泊松比。

重要：所有这些数组中的 N 应该相同。每个观测点恰好有一个三角形和滑移值。

• 输出 disp 是一个形状为 (N, 3) 的数组，包含 x、y、z 方向上的位移分量。
• 输出 strain 是一个形状为 (N, 6) 的数组，表示一个对称张量。 strain[:,0] 是应变 xx 分量，1 是 yy，2 是 zz，3 是 xy，4 是 xz，5 是 yz。

我想得到应力。

使用

stress = cutde.fullspace.strain_to_stress(strain, sm, nu)

从应力转换为应变，假设各向同性线性弹性。 sm 是剪切模量，nu 是泊松比。

所有成对交互矩阵

如果您想创建一个表示每个观测点和每个源三角形的相互作用的矩阵，则有一个不同的接口

from cutde.fullspace import disp_matrix, strain_matrix

disp_mat = disp_matrix(obs_pts, src_tris, nu)
strain_mat = strain_matrix(obs_pts, src_tris, nu)

• obs_pts 是一个形状为 (N_OBS_PTS, 3) 的 np.array。
• src_tris 是一个形状为 (N_SRC_TRIS, 3, 3) 的 np.array，其中第二个维度对应于每个顶点，第三个维度对应于这些顶点的坐标。
• 最后一个参数，nu，是泊松比。
• 输出 disp_mat 是一个形状为 (N_OBS_PTS, 3, N_SRC_TRIS, 3) 的数组。第二个维度对应于观测位移的分量，第四个维度对应于源滑移向量的分量。滑移向量分量的顺序与 disp(...) 和 strain(...) 中的顺序相同。
• 输出 strain_mat 是一个形状为 (N_OBS_PTS, 6, N_SRC_TRIS, 3) 的数组。与上面类似，维度对应于观测应变的分量，顺序与 strain(...) 相同。

注意，要使用此 strain_to_stress 函数与这种矩阵输出，您需要重新排列 strain 数组的轴，以便 6 维应变轴是最后一个轴。您可以使用 np.transpose(...) 来执行此操作。

无矩阵所有成对交互

上述 disp_matrix(...) 产生的矩阵的常见用途是与具有 (N_SRC_TRIS * 3) 个条目和具有 (N_OBS_PTS * 6) 个条目的输入向量进行矩阵-向量乘法。但是，构建整个矩阵可能需要非常大的内存。在某些情况下，在没有计算矩阵本身的情况下计算矩阵-向量乘法可能很有用，这是一种所谓的“无矩阵”操作。为了执行此操作，矩阵项将在需要时重新计算。因此，执行矩阵-向量乘法要慢得多——在我的机器上，大约慢 20 倍。但是，如果内存受限，这种权衡可能值得。

from cutde.fullspace import disp_free, strain_free
disp = disp_free(obs_pts, src_tris, slips, nu)
strain = strain_free(obs_pts, src_tris, slips, nu)

参数与 disp_matrix(...) 相同，另外增加了 slips。 slips 数组是一个形状为 (N_SRC_TRIS, 3) 的数组，包含源滑移向量。

分块交互矩阵

在某些设置中，在计算矩阵的许多子块而不计算整个矩阵时很有用。例如，这对于分层矩阵的近场分量或快速多极近似很有用。

from cutde.fullspace import disp_block, strain_block
disp_matrices, block_idxs = disp_block(
    obs_pts, src_tris, obs_start, obs_end, src_start, src_end, nu
)
strain_matrices, strain_block_idxs = strain_block(
    obs_pts, src_tris, obs_start, obs_end, src_start, src_end, nu
)

• obs_pts、src_tris 和 nu 与 disp_matrix 相同。
• obs_starts 和 obs_end 是具有 N_BLOCKS 个元素的数组，表示每个块中的第一个和最后一个观测点索引。
• src_starts 和 src_end 是具有 N_BLOCKS 个元素的数组，表示每个块中的第一个和最后一个源三角形索引。

输出 disp_matrices 和 strain_matrices 将是密集表示，每个块的边界由 block_idxs 标记。以下是一个提取单个块的示例

disp_matrices, block_idxs = disp_block(obs_pts, src_tris, [0, 5], [5, 10], [0, 2], [2, 4], nu)
block1 = disp_matrices[block_idxs[0]:block_idxs[1]].reshape((5, 3, 2, 3))

自适应交叉逼近（ACA）

有时我们想要计算的矩阵块代表远场交互，观测点都足够远，并且与源三角形分离。在这种情况下，矩阵块近似为低秩。一个近似的矩阵将需要更少的存储空间，并且允许更高效的矩阵向量乘法。自适应交叉逼近是一种用于计算此类低秩表示的算法。参见Grasedyck 2005，以了解对ACA的简单而通用的介绍。或者，参见此处ACA部分，该部分介绍了如何使用逐步使用cutde.fullspace.disp_aca(...)实现。

disp_appxs = cutde.fullspace.disp_aca(
    obs_pts, tris, obs_start, obs_end, src_start, src_end, nu, tol, max_iter
)

与其他所有函数一样，此函数由cutde.fullspace和cutde.halfspace提供。

参数与disp_block(...)相同，增加了tol和max_iter。容差tol以N_BLOCKS长度的数组形式指定，表示真实矩阵与近似矩阵之间的误差矩阵的Frobenius范数。算法不一定达到指定的容差，但应该非常接近。最大迭代次数（等于近似的最大秩）也指定为长度为N_BLOCKS的数组。

输出disp_appxs将是一个表示低秩近似的左和右向量的(U, V)对的列表。为了近似矩阵向量乘积

U, V = disp_appxs[0]
y = U.dot(V.dot(x))

安装

从conda-forge安装更可取，因为应该涉及编译器的问题较少。要从conda-forge安装cutde

conda install -c conda-forge cutde

或者使用pip从PyPI安装

pip install cutde

使用pip安装将通过源代码构建C++扩展，这需要访问非古老版本的GCC或Clang。

这应该足够使用C++/CPU后端。如果您想通过PyCUDA或PyOpenCL使用GPU后端，请按照下面的说明操作。

GPU安装

按照下面的说明安装PyCUDA或PyOpenCL。

PyCUDA

如果您有NVIDIA GPU，请使用以下命令安装PyCUDA

conda config --prepend channels conda-forge
conda install -c conda-forge pycuda

Mac OS X

使用以下命令安装PyOpenCL和PoCL OpenCL驱动程序

conda config --prepend channels conda-forge
conda install pocl pyopencl

Ubuntu + PyOpenCL/PoCL

就像在Mac上一样

conda config --prepend channels conda-forge
conda install pocl pyopencl

Ubuntu + PyOpenCL与系统驱动程序

conda install pyopencl ocl-icd ocl-icd-system

您需要根据您的硬件自行安装系统OpenCL驱动程序。请参阅下面的“其他事项”部分。

Windows

如果您有NVIDIA GPU，请参阅PyCUDA说明。

我不知道是否有人已经在Windows上测试过OpenCL的cutde。它应该不难安装。我预计您将通过conda安装pyopencl，然后安装由您的硬件供应商提供的OpenCL库和驱动程序。请参阅下面的“其他事项”部分。

其他

我建议首先尝试上面与您的系统最相似的说明。如果不起作用，不要担心！几乎在所有最近的硬件和典型的操作系统上都可以安装OpenCL。这些说明可能会有所帮助。这些说明可能会有所帮助。如果您在OpenCL安装方面遇到问题，我将很乐意尝试帮助，但我不能保证有用。

为什么我不能使用Apple CPU OpenCL？

您可能会收到以下消息：cutde不支持Apple CPU OpenCL实现，并且未找到其他平台或设备。请参阅cutde README。

Apple OpenCL实现对Intel CPU的OpenCL标准的支持非常差，并导致许多难以解决的错误。相反，请使用PoCL实现。您可以使用以下命令安装它：conda install -c conda-forge pocl。

开发

要开发cutde，请克隆存储库，并根据environment.yml设置conda环境

conda env create

接下来，对于GPU开发，请按照上面的安装说明安装pycuda或pyopencl。

然后，您应该重新生成基于Mehdi Nikhoo的MATLAB代码的基线测试数据。为此，首先安装octave。在Ubuntu上，只需这样做：

sudo apt-get install octave

然后运行

./tests/setup_test_env

这将运行tests/matlab/gen_test_data.m脚本。

最后，为了检查cutde是否正常工作，运行pytest！

架构

模块总结。

• halfspace.py和fullspace.py - 主入口点。这些是非常薄的包装层，提供用户API。
• coordinators.py - 调用CUDA内核的驱动函数。建议从这里开始！
• geometry.py - 几何辅助函数
• common.cu - MATLAB中主要计算内核的半直接翻译。这些被下面的其他CUDA内核调用。
• pairs.cu - 对TDE进行逐对计算的CUDA内核。
• matrix.cu - 构建矩阵的所有对TDE计算的CUDA内核。
• blocks.cu - 块矩阵计算的CUDA内核。
• free.cu - 矩阵-向量乘法计算的CUDA内核。如果想要构建的矩阵太大而无法在内存中保存，则可以使用。
• aca.cu - 自适应交叉近似实现的CUDA内核。
• backend.py - 在CUDA、OpenCL和C++之间抽象的一层。
• gpu_backend.py - CUDA和OpenCL后端的一些辅助函数
• mako_helpers.py - Mako模板的辅助函数。
• cuda.py - PyCUDA后端。
• opencl.py - PyOpenCL后端。
• cpp.py和cutde.cpp_backend - 结合这两个文件提供了一个可移植层，以便CUDA代码可以实际编译为C++并在CPU上运行，尽管速度较慢。

tests/tde_profile.py脚本用于评估性能。

一些测试被标记为慢。要运行这些测试，运行pytest --runslow。

如果您有多个后端可用并安装了cutde，则将首选CUDA，然后是OpenCL，最后回退到C++后端。如果您希望指定要使用哪个后端，可以将环境变量CUTDE_USE_BACKEND设置为cuda、opencl或cpp。

README.md是从docs/中的模板自动生成的。要运行此过程，请运行docs/build_docs。