feat: add GCD parallel implementation for macOS

+utils/computeFourthOrderDiff.m

@@ -0,0 +1,90 @@
+%% 四阶中心差分计算空间导数
+% 功能：计算二维场的空间导数，内部使用四阶中心差分，边界使用二阶差分
+% 
+% 说明：
+%   1. 主要功能：
+%      - 计算x方向空间导数
+%      - 计算y方向空间导数
+%      - 边界特殊处理
+%   2. 计算方法：
+%      - 内部点：四阶中心差分
+%      - 边界点：二阶前向/后向差分
+%      - 次边界点：二阶中心差分
+%
+% 使用方法：
+%   [dx, dy] = utils.computeFourthOrderDiff(field, deltax, deltay)
+%
+% 输入参数：
+%   - field: 输入场 (2D矩阵)
+%   - deltax: x方向网格间距
+%   - deltay: y方向网格间距
+%
+% 输出参数：
+%   - dx: x方向导数
+%   - dy: y方向导数
+%
+% 依赖项：
+%   - 无
+%
+% 注意事项：
+%   - 输入场必须是二维矩阵
+%   - 网格间距必须为正数
+%   - 边界处理使用较低阶差分以保持稳定性
+%
+% 作者：StarrMoonn
+% 日期：2025-01-16
+%
+
+function [dx, dy] = computeFourthOrderDiff(field, deltax, deltay)
+    % 四阶中心差分计算空间导数
+    %
+    % 输入参数:
+    %   field: 输入场 (2D矩阵)
+    %   deltax: x方向网格间距
+    %   deltay: y方向网格间距
+    %
+    % 输出参数:
+    %   dx: x方向导数
+    %   dy: y方向导数
+    %
+    % 说明：
+    %   - 内部点使用四阶中心差分
+    %   - 边界点使用二阶差分
+    %   - 适用于地震波场正演和反演计算
+
+    % 获取场的大小
+    [nx, ny] = size(field);
+
+    % 初始化导数数组
+    dx = zeros(nx, ny);
+    dy = zeros(nx, ny);
+
+    % 四阶中心差分系数（修正系数）
+    c1 = -8/12;   % 相邻点系数
+    c2 = 1/12;    % 间隔一个点的系数
+
+    % 计算x方向导数（四阶中心差分）
+    for i = 3:nx-2
+        dx(i,:) = (c2*(field(i-2,:) - field(i+2,:)) + ...
+                   c1*(field(i-1,:) - field(i+1,:))) / deltax;
+    end
+
+    % 计算y方向导数（四阶中心差分）
+    for j = 3:ny-2
+        dy(:,j) = (c2*(field(:,j-2) - field(:,j+2)) + ...
+                   c1*(field(:,j-1) - field(:,j+1))) / deltay;
+    end
+
+    % 处理边界（使用二阶差分）
+    % x方向边界
+    dx(1,:) = (-3*field(1,:) + 4*field(2,:) - field(3,:)) / (2*deltax);  % 二阶前向差分
+    dx(2,:) = (field(3,:) - field(1,:)) / (2*deltax);  % 二阶中心差分
+    dx(nx-1,:) = (field(nx,:) - field(nx-2,:)) / (2*deltax);  % 二阶中心差分
+    dx(nx,:) = (field(nx-2,:) - 4*field(nx-1,:) + 3*field(nx,:)) / (2*deltax);  % 二阶后向差分
+
+    % y方向边界
+    dy(:,1) = (-3*field(:,1) + 4*field(:,2) - field(:,3)) / (2*deltay);  % 二阶前向差分
+    dy(:,2) = (field(:,3) - field(:,1)) / (2*deltay);  % 二阶中心差分
+    dy(:,ny-1) = (field(:,ny) - field(:,ny-2)) / (2*deltay);  % 二阶中心差分
+    dy(:,ny) = (field(:,ny-2) - 4*field(:,ny-1) + 3*field(:,ny)) / (2*deltay);  % 二阶后向差分
+end 

+utils/computeGradientField.m

@@ -0,0 +1,29 @@
+%% MATLAB的gradient函数计算场的空间导数
+%
+% 功能：
+%   二阶中心差分计算空间导数，使用MATLAB内置的gradient函数
+%
+% 输入参数：
+%   field: 输入场 (2D矩阵)
+%   deltax: x方向网格间距 (单位：m)
+%   deltay: y方向网格间距 (单位：m)
+%
+% 输出参数：
+%   dx: x方向导数 (与输入场同维度)
+%   dy: y方向导数 (与输入场同维度)
+%
+% 注意事项：
+%   - gradient函数返回顺序为[dy, dx]，此函数将其调整为[dx, dy]
+%   - 导数已除以网格间距，具有正确的物理单位
+%
+% 作者：StarrMoonn
+% 日期：2025-01-16
+%
+function [dx, dy] = computeGradientField(field, deltax, deltay)
+    % 计算导数并除以实际间距
+    [dy, dx] = gradient(field);
+
+    % 除以网格间距得到物理单位
+    dx = dx / deltax;
+    dy = dy / deltay;
+end 

+utils/compute_wave_propagation_cpu1.m

@@ -0,0 +1,98 @@
+function [vx, vy] = compute_wave_propagation_cpu1(obj)
+    % 从对象中获取计算参数
+    dx = obj.DELTAX;
+    dy = obj.DELTAY;
+    dt = obj.DELTAT;
+
+    % 获取当前状态变量的引用（为了代码简洁性）
+    vx = obj.vx;
+    vy = obj.vy;
+    sigmaxx = obj.sigmaxx;
+    sigmayy = obj.sigmayy;
+    sigmaxy = obj.sigmaxy;
+
+    % 计算应力场 sigmaxx 和 sigmayy
+    for j = 2:obj.NY
+        for i = 1:obj.NX-1
+            % 计算速度梯度
+            value_dvx_dx = (vx(i+1,j) - vx(i,j)) / dx;
+            value_dvy_dy = (vy(i,j) - vy(i,j-1)) / dy;
+
+            % 更新PML记忆变量
+            obj.memory_dvx_dx(i,j) = obj.b_x_half(i) * obj.memory_dvx_dx(i,j) + ...
+                                    obj.a_x_half(i) * value_dvx_dx;
+            obj.memory_dvy_dy(i,j) = obj.b_y(j) * obj.memory_dvy_dy(i,j) + ...
+                                    obj.a_y(j) * value_dvy_dy;
+
+            % 计算最终值
+            value_dvx_dx = value_dvx_dx / obj.K_x_half(i) + obj.memory_dvx_dx(i,j);
+            value_dvy_dy = value_dvy_dy / obj.K_y(j) + obj.memory_dvy_dy(i,j);
+
+            % 更新应力
+            sigmaxx(i,j) = sigmaxx(i,j) + dt * (...
+                obj.c11(i,j) * value_dvx_dx + ...
+                obj.c13(i,j) * value_dvy_dy);
+
+            sigmayy(i,j) = sigmayy(i,j) + dt * (...
+                obj.c13(i,j) * value_dvx_dx + ...
+                obj.c33(i,j) * value_dvy_dy);
+        end
+    end
+
+    % 计算剪应力 sigmaxy
+    for j = 1:obj.NY-1
+        for i = 2:obj.NX
+            value_dvy_dx = (vy(i,j) - vy(i-1,j)) / dx;
+            value_dvx_dy = (vx(i,j+1) - vx(i,j)) / dy;
+
+            obj.memory_dvy_dx(i,j) = obj.b_x(i) * obj.memory_dvy_dx(i,j) + ...
+                                    obj.a_x(i) * value_dvy_dx;
+            obj.memory_dvx_dy(i,j) = obj.b_y_half(j) * obj.memory_dvx_dy(i,j) + ...
+                                    obj.a_y_half(j) * value_dvx_dy;
+
+            value_dvy_dx = value_dvy_dx / obj.K_x(i) + obj.memory_dvy_dx(i,j);
+            value_dvx_dy = value_dvx_dy / obj.K_y_half(j) + obj.memory_dvx_dy(i,j);
+
+            sigmaxy(i,j) = sigmaxy(i,j) + ...
+                obj.c44(i,j) * (value_dvy_dx + value_dvx_dy) * dt;
+        end
+    end
+
+    % 计算x方向速度场
+    for j = 2:obj.NY
+        for i = 2:obj.NX
+            value_dsigmaxx_dx = (sigmaxx(i,j) - sigmaxx(i-1,j)) / dx;
+            value_dsigmaxy_dy = (sigmaxy(i,j) - sigmaxy(i,j-1)) / dy;
+
+            obj.memory_dsigmaxx_dx(i,j) = obj.b_x(i) * obj.memory_dsigmaxx_dx(i,j) + ...
+                                         obj.a_x(i) * value_dsigmaxx_dx;
+            obj.memory_dsigmaxy_dy(i,j) = obj.b_y(j) * obj.memory_dsigmaxy_dy(i,j) + ...
+                                         obj.a_y(j) * value_dsigmaxy_dy;
+
+            value_dsigmaxx_dx = value_dsigmaxx_dx / obj.K_x(i) + obj.memory_dsigmaxx_dx(i,j);
+            value_dsigmaxy_dy = value_dsigmaxy_dy / obj.K_y(j) + obj.memory_dsigmaxy_dy(i,j);
+
+            vx(i,j) = vx(i,j) + ...
+                (value_dsigmaxx_dx + value_dsigmaxy_dy) * dt / obj.rho(i,j);
+        end
+    end
+
+    % 计算y方向速度场
+    for j = 1:obj.NY-1
+        for i = 1:obj.NX-1
+            value_dsigmaxy_dx = (sigmaxy(i+1,j) - sigmaxy(i,j)) / dx;
+            value_dsigmayy_dy = (sigmayy(i,j+1) - sigmayy(i,j)) / dy;
+
+            obj.memory_dsigmaxy_dx(i,j) = obj.b_x_half(i) * obj.memory_dsigmaxy_dx(i,j) + ...
+                                         obj.a_x_half(i) * value_dsigmaxy_dx;
+            obj.memory_dsigmayy_dy(i,j) = obj.b_y_half(j) * obj.memory_dsigmayy_dy(i,j) + ...
+                                         obj.a_y_half(j) * value_dsigmayy_dy;
+
+            value_dsigmaxy_dx = value_dsigmaxy_dx / obj.K_x_half(i) + obj.memory_dsigmaxy_dx(i,j);
+            value_dsigmayy_dy = value_dsigmayy_dy / obj.K_y_half(j) + obj.memory_dsigmayy_dy(i,j);
+
+            vy(i,j) = vy(i,j) + ...
+                (value_dsigmaxy_dx + value_dsigmayy_dy) * dt / obj.rho(i,j);
+        end
+    end
+end 

+utils/compute_wave_propagation_cpu2.m

@@ -0,0 +1,59 @@
+% +utils/compute_wave_propagation.m
+function [vx, vy] = compute_wave_propagation_cpu2(obj)
+    % 从对象中获取所需参数
+    vx = obj.vx;           % x方向速度场
+    vy = obj.vy;           % y方向速度场
+    sigmaxx = obj.sigmaxx; % x方向应力
+    sigmayy = obj.sigmayy; % y方向应力
+    sigmaxy = obj.sigmaxy; % 剪切应力
+
+    % 内存变量
+    memory_dvx_dx = obj.memory_dvx_dx;
+    memory_dvy_dy = obj.memory_dvy_dy;
+    memory_dvy_dx = obj.memory_dvy_dx;
+    memory_dvx_dy = obj.memory_dvx_dy;
+    memory_dsigmaxx_dx = obj.memory_dsigmaxx_dx;
+    memory_dsigmaxy_dy = obj.memory_dsigmaxy_dy;
+    memory_dsigmaxy_dx = obj.memory_dsigmaxy_dx;
+    memory_dsigmayy_dy = obj.memory_dsigmayy_dy;
+
+    % 材料参数
+    c11 = obj.c11;
+    c13 = obj.c13;
+    c33 = obj.c33;
+    c44 = obj.c44;
+    rho = obj.rho;
+
+    % PML参数
+    b_x = obj.b_x;
+    b_y = obj.b_y;
+    b_x_half = obj.b_x_half;
+    b_y_half = obj.b_y_half;
+    a_x = obj.a_x;
+    a_y = obj.a_y;
+    a_x_half = obj.a_x_half;
+    a_y_half = obj.a_y_half;
+    K_x = obj.K_x;
+    K_y = obj.K_y;
+    K_x_half = obj.K_x_half;
+    K_y_half = obj.K_y_half;
+
+    % 计算参数
+    DELTAX = obj.DELTAX;
+    DELTAY = obj.DELTAY;
+    DELTAT = obj.DELTAT;
+    NX = obj.NX;
+    NY = obj.NY;
+
+    % 调用编译好的 MEX 文件
+    [vx, vy] = compute_wave_propagation(vx, vy, sigmaxx, sigmayy, sigmaxy, ...
+        memory_dvx_dx, memory_dvy_dy, memory_dvy_dx, memory_dvx_dy, ...
+        memory_dsigmaxx_dx, memory_dsigmaxy_dy, memory_dsigmaxy_dx, memory_dsigmayy_dy, ...
+        c11, c13, c33, c44, rho, ...
+        b_x, b_y, b_x_half, b_y_half, ...
+        a_x, a_y, a_x_half, a_y_half, ...
+        K_x, K_y, K_x_half, K_y_half, ...
+        DELTAX, DELTAY, DELTAT, NX, NY);
+end
+
+

+utils/compute_wave_propagation_cpu3.m

@@ -0,0 +1,56 @@
+function [vx, vy] = compute_wave_propagation_cpu3(obj)
+    % 从对象中获取所需参数
+    vx = obj.vx;           % x方向速度场
+    vy = obj.vy;           % y方向速度场
+    sigmaxx = obj.sigmaxx; % x方向应力
+    sigmayy = obj.sigmayy; % y方向应力
+    sigmaxy = obj.sigmaxy; % 剪切应力
+
+    % 内存变量
+    memory_dvx_dx = obj.memory_dvx_dx;
+    memory_dvy_dy = obj.memory_dvy_dy;
+    memory_dvy_dx = obj.memory_dvy_dx;
+    memory_dvx_dy = obj.memory_dvx_dy;
+    memory_dsigmaxx_dx = obj.memory_dsigmaxx_dx;
+    memory_dsigmaxy_dy = obj.memory_dsigmaxy_dy;
+    memory_dsigmaxy_dx = obj.memory_dsigmaxy_dx;
+    memory_dsigmayy_dy = obj.memory_dsigmayy_dy;
+
+    % 材料参数
+    c11 = obj.c11;
+    c13 = obj.c13;
+    c33 = obj.c33;
+    c44 = obj.c44;
+    rho = obj.rho;
+
+    % PML参数
+    b_x = obj.b_x;
+    b_y = obj.b_y;
+    b_x_half = obj.b_x_half;
+    b_y_half = obj.b_y_half;
+    a_x = obj.a_x;
+    a_y = obj.a_y;
+    a_x_half = obj.a_x_half;
+    a_y_half = obj.a_y_half;
+    K_x = obj.K_x;
+    K_y = obj.K_y;
+    K_x_half = obj.K_x_half;
+    K_y_half = obj.K_y_half;
+
+    % 计算参数
+    DELTAX = obj.DELTAX;
+    DELTAY = obj.DELTAY;
+    DELTAT = obj.DELTAT;
+    NX = obj.NX;
+    NY = obj.NY;
+
+    % 调用编译好的 OpenMP MEX 文件
+    [vx, vy] = compute_wave_propagation_gcd(vx, vy, sigmaxx, sigmayy, sigmaxy, ...
+        memory_dvx_dx, memory_dvy_dy, memory_dvy_dx, memory_dvx_dy, ...
+        memory_dsigmaxx_dx, memory_dsigmaxy_dy, memory_dsigmaxy_dx, memory_dsigmayy_dy, ...
+        c11, c13, c33, c44, rho, ...
+        b_x, b_y, b_x_half, b_y_half, ...
+        a_x, a_y, a_x_half, a_y_half, ...
+        K_x, K_y, K_x_half, K_y_half, ...
+        DELTAX, DELTAY, DELTAT, NX, NY);
+end 

README.md

@@ -1,6 +1,6 @@
 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;蝶恋花<br>
 
-&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;阅尽天涯离别苦， 不道归来，零落花如许。花底相看无一语， 绿窗春与天俱暮。<br>
+&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;阅尽天涯离别苦， 不道归来，零落花如许。花底相看无一语， 绿窗春与天俱莫。<br>
 
 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;待把相思灯下诉， 一缕新欢，旧恨千千缕。最是人间留不住， 朱颜辞镜花辞树。