障碍似乎只在一个时间窗口内同步

Posted 2023-03-17

【中文标题】障碍似乎只在一个时间窗口内同步

【英文标题】：Barriers seems to synchronize only within a time window

【发布时间】：2019-08-03 13:55:15

【问题描述】：

目前我尝试实施 FDTD 方法来求解麦克斯韦方程组 使用 OpenCL。算法很简单，计算当前的h-field 从旧电场，并从当前的 h 场计算当前的 e 场。然后开始下一次迭代。在 OpenCL 中出现同步问题 因为在计算 h 场之前无法计算 e 场，并且必须同步开始下一次迭代。因此只需要使用一个工作组。我通过使全局和局部工作空间大小相同来确保这一点。 我的 gpu 的最大工作项值为 256，这意味着对于 3d 问题空间，我可以在每个维度上拥有 6 个工作项。这就是为什么每个工作项都计算一个字段值的立方体。我遇到的问题是，只要每个工作项在每个维度中计算的字段值少于 12 个，总共 1728 个字段值，算法就会产生正确的结果。 如果工作项必须计算更多的字段值，则算法的结果会越来越差。我也用双精度浮点测试了算法，结果更糟。 我对双精度结果的怀疑是，所有字段值的计算都需要很多时间，并且同步不再正常工作。以下是 OpenCL 内核的源代码：

内核无效 fdtd3d_noiter_fp32( 全局 float3* old_e, 全局 float3* old_h, 全局 float3* current_e, 全局 float3* current_h, 全局 float3* e_curl_integral, 全局 float3* e_field_integral, 全局 float3* h_curl_integral, 全局 float3* h_field_integral, 常量 float4* ex_factor, 常数 float4* ey_factor, 常量 float4* ez_factor, 常量 float4* hx_factor, 常量 float4* hy_factor, 常量 float4* hz_factor, 常量 char* 几何， 浮动 grid_width_x, 浮动grid_width_y， 浮动grid_width_z， 整数宽度， 深入了解， 浮动最大时间， 浮动时间步长， int batch_size_x, int batch_size_y, int batch_size_z, 常量 char* source_factor, 浮动源振幅， 浮动 source_ramp_len， 浮动源频率） // 获取实际id 常量 int ID = get_global_id(0) * batch_size_x + get_global_id(1) * batch_size_y * uwidth + get_global_id(2) * batch_size_z * uwidth * udepth + 1 + uwidth + uwidth * udepth； int id = ID; for (float current_time = 0; current_time < max_time; current_time += time_step) for (int z = 0; z < batch_size_z; z++) for (int y = 0; y < batch_size_y; y++) for (int x = 0; x < batch_size_x; x++) id = ID + x + y * uwidth + z * uwidth * udepth; // calculate the current magnetic field from the old electric field float hx_curl_term = ((old_e[id + uwidth].z - old_e[id].z) / grid_width_y) - ((old_e[id + uwidth * udepth].y - old_e[id].y) / grid_width_z); float hy_curl_term = ((old_e[id + uwidth * udepth].x - old_e[id].x) / grid_width_z) - ((old_e[id + 1].z - old_e[id].z) / grid_width_x); float hz_curl_term = ((old_e[id + 1].y - old_e[id].y) / grid_width_x) - ((old_e[id + uwidth].x - old_e[id].x) / grid_width_y); h_curl_integral[id].x += hx_curl_term; h_curl_integral[id].y += hy_curl_term; h_curl_integral[id].z += hz_curl_term; h_field_integral[id] += old_h[id]; current_h[id].x = + hx_factor[id].x * old_h[id].x - hx_factor[id].y * hx_curl_term - hx_factor[id].z * h_curl_integral[id].x - hx_factor[id].w * h_field_integral[id].x; current_h[id].y = + hy_factor[id].x * old_h[id].y - hy_factor[id].y * hy_curl_term - hy_factor[id].z * h_curl_integral[id].y - hy_factor[id].w * h_field_integral[id].y; current_h[id].z = + hz_factor[id].x * old_h[id].z - hz_factor[id].y * hz_curl_term - hz_factor[id].z * h_curl_integral[id].z - hz_factor[id].w * h_field_integral[id].z; barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE); // calculate the current electric field from the current magnetic field for (int z = 0; z < batch_size_z; z++) for (int y = 0; y < batch_size_y; y++) for (int x = 0; x < batch_size_x; x++) id = ID + x + y * uwidth + z * uwidth * udepth; float ex_curl_term = ((current_h[id].z - current_h[id - uwidth].z) / grid_width_y) - ((current_h[id].y - current_h[id - uwidth * udepth].y) / grid_width_z); float ey_curl_term = ((current_h[id].x - current_h[id - uwidth * udepth].x) / grid_width_z) - ((current_h[id].z - current_h[id - 1].z) / grid_width_x); float ez_curl_term = ((current_h[id].y - current_h[id - 1].y) / grid_width_x) - ((current_h[id].x - current_h[id - uwidth].x) / grid_width_y); e_curl_integral[id].x += ex_curl_term; e_curl_integral[id].y += ey_curl_term; e_curl_integral[id].z += ez_curl_term; e_field_integral[id] += old_e[id]; current_e[id].x = ex_factor[id].x * old_e[id].x - ex_factor[id].y * e_field_integral[id].x + ex_factor[id].z * ex_curl_term + ex_factor[id].w * e_curl_integral[id].x; current_e[id].y = ey_factor[id].x * old_e[id].y - ey_factor[id].y * e_field_integral[id].y + ey_factor[id].z * ey_curl_term + ey_factor[id].w * e_curl_integral[id].y; current_e[id].z = ez_factor[id].x * old_e[id].z - ez_factor[id].y * e_field_integral[id].z + ez_factor[id].z * ez_curl_term + ez_factor[id].w * e_curl_integral[id].z; // zero ez out if the soure_factor is 1 current_e[id].z -= abs(source_factor[id]) * current_e[id].z; // than add the source. this way a hard source is created! current_e[id].z += source_factor[id] * source_amplitude * ramped_sin_fp32(current_time, source_ramp_len, source_frequency); current_e[id] *= geometry[id]; // synchronize the workitems after calculating the current fields barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE); // now memory consistency is assured and the fields can be copied for (int z = 0; z < batch_size_z; z++) for (int y = 0; y < batch_size_y; y++) for (int x = 0; x < batch_size_x; x++) id = ID + x + y * uwidth + z * uwidth * udepth; old_h[id] = current_h[id]; old_e[id] = current_e[id]; barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);

我希望代码不是难以理解，有人可以帮助我！

躺着

【参考方案1】：

让我们从小事开始：

屏障（CLK_LOCAL_MEM_FENCE

如果不使用局部变量，则使用局部屏障毫无意义。

我对双精度结果的怀疑是，所有字段值的计算需要很长时间，并且同步不再正常工作

GPU 中的同步是由硬件完成的，虽然可能，但问题不大。

现在是个大问题：

在 OpenCL 中出现同步问题，因为 e-field 不能 在计算 h 场和下一次迭代之前计算 必须同步启动。因此只需一个工作组 用过的。我通过使全局和局部工作空间大小相同来确保这一点。

...这是绝对不行的。您实际上是在序列化您的代码。这会降低性能并导致 GPU 利用率非常低，而且运行速度可能比在 CPU 上慢。

我没有仔细查看您的代码，但有可能 (AFAICT) 您不必使用单个工作组。我写它的方式（伪代码）：

OpenCL 代码：

kernel1(arguments)  calculate the current magnetic field from the old electric field 
kernel2(arguments)  calculate the current electric field from the current magnetic field 
kernel3(arguments)  copy current to old (both e/h) 

主机代码：

for (float current_time = 0; current_time < max_time; current_time += time_step)
   
       clEnqueueNDRangeKernel(queue, kernel1, ...);
       clEnqueueNDRangeKernel(queue, kernel2, ...);
       clEnqueueNDRangeKernel(queue, kernel3, ...);
   

这是否是一个好主意，取决于很多事情 - 问题大小、时间步长等。但它允许您使用 3D 网格进行全局工作，并且 GPU 将并行运行这些。也许尝试在 github 上找到一些计算 Maxwell 方程的 OpenCL 或 CUDA 代码。

【讨论】：

感谢您详细更正我的代码，即使我仍然好奇导致这种行为的原因，我也会按照您的建议重写我的代码。