vlsplt 串并行 Van Leer 差异对照

输入范围

dyn3d\vlsplt.F: L1-1135（串行主例程 + vlx/vly/vlz + minmaxq）
dyn3dpar\vlsplt_p.F: L1-1340（并行主例程 + vlx_p/vly_p/vlz_p + minmaxq_p）

Mars 运行参与度：条件经过。Mars 串行构建经过 vlsplt，并行构建经过 vlsplt_p。两者数值框架相同但执行模型和数据布局有系统性差异。本页面向需要同时理解两套代码的维护者和复现者。

页面定位

本页逐项对照 vlsplt（串行）和 vlsplt_p（并行）的结构差异、数值等价性和潜在分歧点。各方向的详细算法说明见 vlsplt-serial-split-kernels 和 vlsplt-p-halo-split。

架构差异总览

维度	串行 vlsplt	并行 vlsplt_p
计算范围	全域 `iip2:ip1jm`	band `ij_begin:ij_end`
域分解	无	`bands.F90` SetDistrib
通信	无	Register_Hallo + SendRequest/WaitRecvRequest/WaitSendRequest
通信重叠	无	边界带先算 → halo 发送 → 内区计算 → 等待
线程	无 OpenMP	OMP DO SCHEDULE(STATIC,OMP_CHUNK)
线程私有	无 THREADPRIVATE	vly_p 几何量和计时变量 + vlz_p 的 `first/testcpu/temps*`
工作数组	栈上固定维度	vlz_p 的 wq/Ratio 动态 ALLOCATABLE SAVE
split 调用次数	每方向 1 次全域	vlx/vlz 各 3 次分段（边界+内区），vly 2 次全域
源文件位置	`libf/dyn3d/`	`libf/dyn3dpar/`

主例程对照

签名

! 串行
SUBROUTINE vlsplt(q, pente_max, masse, u_m, w, iq)

! 并行
SUBROUTINE vlsplt_p(q, pente_max, masse, u_m, w, iq)

签名完全一致。差异全在内部实现。

垂直通量维度

! 串行
REAL mw(ip1jmp1, llm+1)

! 并行
REAL mw(ip1jmp1, llm+1, nqtot)

并行版为每个 tracer 保留独立的垂直通量 mw(:,:,iq)。这是因为并行主例程对每个 iq 分别调用五步 split，而不是在主例程内循环 iq——后者在串行版中允许复用同一 mw 数组。

五步 split 顺序

两个版本的 Strang split 顺序相同：x → y → z → y → x。但执行方式不同：

串行：vlx(全域) → vly(全域) → vlz(全域) → vly(全域) → vlx(全域)

并行：
  vlx_p(左边界) + vlx_p(右边界) → Register_Hallo → SendRequest → vlx_p(内区) → WaitRecv
  vly_p(全域 band)
  vlz_p(左边界) + vlz_p(右边界) → Register_Hallo → SendRequest → vlz_p(内区) → WaitRecv
  vly_p(全域 band)
  vlx_p(全域 band)
  WaitSendRequest

关键差异：并行版的 vlx 和 vlz 被拆成"边界带 + halo + 内区"三段以重叠计算和通信；vly 不需要 halo 因为 N-S 数据依赖已在 band 内部。

通量缩放

! 串行：每步前缩放
u_m = u_m * 0.5
w   = w   * 0.5
! 最后一步 vlx 前恢复
u_m = u_m * 2.

! 并行：相同逻辑
u_m = u_m * 0.5
w   = w   * 0.5
! 最后一步前恢复
u_m = u_m * 2.

通量缩放策略完全相同。两个版本的中间步使用半通量、最后步恢复全通量。

q 回写

! 串行
DO iq = 1, nqtot
   q(:,:,iq) = zq(:,:)
ENDDO

! 并行
! 每个 iq 的五步 split 结束后直接写回
q(ij,l,iq) = zq(ij,l)

串行版在主例程内用临时 zq 做完所有 tracer 后统一回写；并行版在每个 tracer 的五步 split 完成后立即回写。

vlx 对照

签名

! 串行
RECURSIVE SUBROUTINE vlx(q, pente_max, masse, u_m, iq)

! 并行
RECURSIVE SUBROUTINE vlx_p(q, pente_max, masse, u_m, ijb_x, ije_x, iq)

并行版增加 ijb_x/ije_x 两个参数，使主例程可以分段调用（边界带或内区）。串行版硬编码全域 iip2:ip1jm。

极点行处理

! 串行：无显式排除（极点行在全域内自然包含）

! 并行：
if (pole_nord .and. ijb==1) ijb = ijb + iip1
if (pole_sud .and. ije==ip1jmp1) ije = ije - iip1

并行版在拥有极点的进程上跳过极点行，因为极点没有确定的经度方向。串行版不需要这个判断因为全域计算隐含了极点行的正确索引。

数值计算

Van Leer 斜率限制、上游通量公式和 Semi-Lagrangian CFL > 1 累积逻辑在两个版本中完全相同。包括：

标准斜率：pente_max > 0 时使用 0.5*(dxq(i-1)+dxq(i))，受 pente_max*min(|dxq(i-1)|,|dxq(i)|) 限制。
乘积斜率：pente_max < -1e-5 时使用 sign(min(|p1|,|p2|),p1) 其中 p1=2*dxq(i-1)*dxq(i)/(dxq(i-1)+dxq(i))。
Semi-Lagrangian 累积：zdum < 0 时逐格点追踪通量累积。

E-W 周期包裹

两个版本使用相同的 mod(i-2+iim,iim)+1 公式做 E-W 周期包裹。并行版在 band 边界处不需要包裹（halo 提供邻居数据），但在 band 内部和极点附近的周期行仍然使用相同公式。

父子 tracer

递归 vlx/vlx_p 的父子比例输送逻辑相同：masseq = max(masse*q, masseqmin)、Ratio = q_child/q_parent（受 qperemin 阈值保护）、递归调用、q_child = q_parent * Ratio 恢复。并行版的所有循环限定在 [ijb, ije] 范围内。

vly 对照

签名

! 串行
RECURSIVE SUBROUTINE vly(q, pente_max, masse, masse_adv_v, iq)

! 并行
RECURSIVE SUBROUTINE vly_p(q, pente_max, masse, masse_adv_v, iq)

签名一致。但并行版不接收 band 参数，而是直接使用 parallel_lmdz 的 ij_begin/ij_end。

几何缓存变量

! 串行：
SAVE first, testcpu, temps0-5
! 隐式 SAVE（DATA 初始化）：sinlon, coslon, sinlondlon, coslondlon, airej2, airejjm

! 并行：
SAVE first, testcpu, temps0-5
!$OMP THREADPRIVATE(first, testcpu, temps0-5)
SAVE sinlon, coslon, sinlondlon, coslondlon
!$OMP THREADPRIVATE(sinlon, coslon, sinlondlon, coslondlon)
SAVE airej2, airejjm
!$OMP THREADPRIVATE(airej2, airejjm)

并行版需要显式 THREADPRIVATE 声明，因为 OpenMP 的默认共享语义会让这些 SAVE 变量成为所有线程共享的单一副本。串行版不需要 THREADPRIVATE 因为没有 OpenMP。

扩展 band 范围

! 串行：无扩展，使用全域
DO ij = iip2, ip1jm   ! v 点斜率
DO ij = iip2, ip1jm   ! 标量点斜率

! 并行：需要扩展以获取 N-S 邻居
ijb = ij_begin - 2*iip1    ! v 点斜率
ije = ij_end + iip1
ijb = ij_begin - iip1      ! 标量点斜率
ije = ij_end + iip1

并行版的 band 扩展数据来自之前 vlx_p 阶段的 halo 交换结果。串行版不需要扩展因为整个域都在内存中。

极点处理

这是串并行差异最显著的部分：

! 串行（vly L586-606）：
! 无条件计算两个极点
qpns = SUM(...) / SUM(...)     ! 北极平均
qpsn = SUM(...) / SUM(...)     ! 南极平均
! sin/cos 滤波
GOTO 8888                      ! 无条件跳过限制器
8888 CONTINUE
DO ij = 1, iip1                ! 极点 dyq 归零
   dyq(ij,l) = 0.
   dyq(ip1jm+ij,l) = 0.
ENDDO

! 并行（vly_p）：
if (pole_nord) then             ! 仅北极进程
   qpns = SUM(...) / SUM(...)
   ! sin/cos 滤波（相同逻辑）
   DO ij = 1, iip1
      dyq(ij,l) = 0.
   ENDDO
endif
if (pole_sud) then              ! 仅南极进程
   qpsn = SUM(...) / SUM(...)
   ! sin/cos 滤波（相同逻辑）
   DO ij = ip1jm+1, ip1jmp1
      dyq(ij,l) = 0.
   ENDDO
endif

关键差异：

串行版无条件处理两个极点；并行版按 pole_nord/pole_sud 条件处理，每个进程只处理自己拥有的极点。
串行版使用 GOTO 8888 无条件跳过 fn/fs 比例限制器段（L662-L676）后归零极点 dyq。并行版注释掉了 GOTO 8888，直接在内条件块中归零。
两者最终效果相同：极点 dyq 被归零，极点 Van Leer 修正被禁用。
注释 c ym tout cela ne sert pas a grand chose（vly_p）确认 sin/cos 投影计算是残留代码，因为 dyq 最终被归零。

极点 SSUM 累积

! 串行：无条件
convpn  = SSUM(iim, qbyv(1,l), 1)
convmpn = SSUM(iim, masse_adv_v(1,l), 1)
! ... 北极更新 ...
convps  = -SSUM(iim, qbyv(ip1jm-iim,l), 1)
! ... 南极更新 ...

! 并行：条件化
if (pole_nord) then
   convpn  = SSUM(iim, qbyv(1,l), 1)
   ! ... 北极更新（相同公式）...
endif
if (pole_sud) then
   convps  = -SSUM(iim, qbyv(ip1jm-iim,l), 1)
   ! ... 南极更新（相同公式）...
endif

公式完全相同，只是加了进程位置条件。

vlz 对照

签名

! 串行
RECURSIVE SUBROUTINE vlz(q, pente_max, masse, w, iq)

! 并行
RECURSIVE SUBROUTINE vlz_p(q, pente_max, masse, w, ijb_x, ije_x, iq)

并行版增加 ijb_x/ije_x 参数。

工作数组

! 串行：
REAL wq(ip1jmp1, llm+1)
REAL dzq(ip1jmp1, llm), dzqw(ip1jmp1, llm), adzqw(ip1jmp1, llm)

! 并行：
REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: wq     ! (ip1jmp1, llm+1, nqtot)
REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: Ratio   ! (ip1jmp1, llm, nqtot)
REAL, SAVE :: dzq(ip1jmp1, llm), dzqw(ip1jmp1, llm), adzqw(ip1jmp1, llm)

主要差异：

wq 从二维 (ip1jmp1,llm+1) 扩展为三维 (ip1jmp1,llm+1,nqtot)，每个 tracer 保留独立通量。
新增 Ratio 三维数组，用于父子 tracer 比例恢复。
使用 ALLOCATABLE, SAVE 而非栈上固定数组，支持延迟分配。

分配模式

! 串行：直接声明，编译器在栈或 BSS 段分配

! 并行：
IF (first) THEN
   first = .FALSE.
   !$OMP MASTER
   ALLOCATE(wq(ip1jmp1, llm+1, nqtot))
   ALLOCATE(Ratio(ip1jmp1, llm, nqtot))
   !$OMP END MASTER
   !$OMP BARRIER
END IF

并行版使用 OMP MASTER + BARRIER 确保只有主线程分配，所有线程等待完成后才使用数组。注释 ! vlz_loc si on veut qu'elles soient vues par tous les threads 解释了为什么不用 PRIVATE。

垂直边界条件

! 串行：
DO ij = 1, ip1jmp1
   dzq(ij,1) = 0.
   dzq(ij,llm) = 0.
   wq(ij,llm+1) = 0.
   wq(ij,1) = 0.
ENDDO

! 并行：
!$OMP MASTER
DO ij = ijb, ije
   dzq(ij,1) = 0.
   dzq(ij,llm) = 0.
ENDDO
!$OMP END MASTER
!$OMP BARRIER
! ...（通量计算后）...
!$OMP MASTER
DO ij = ijb, ije
   wq(ij,llm+1,iq) = 0.
   wq(ij,1,iq) = 0.
ENDDO
!$OMP END MASTER
!$OMP BARRIER

串行版设置全域边界，并行版只设置 band 范围内 [ijb, ije] 的边界。并行版用 OMP MASTER 确保只有一个线程设置，BARRIER 确保其他线程等待。

父子 tracer 通量传递

! 串行（vlz L976-981）：
masse(ij,l,iq2) = masseq
Ratio(ij,l,iq2) = q_child / q_parent
wq(ij,l,iq2) = wq(ij,l,iq)    ! 注意：子代 wq 直接复制父代 wq

! 并行（vlz_p L1209-1216）：
masse(ij,l,iq2) = masseq
Ratio(ij,l,iq2) = q_child / q_parent
w(ij,l,iq2) = wq(ij,l,iq)     ! 注意：写入 w 而非 wq
! 注释：correction bug le 15mai2015

重要差异：串行版将父代 wq 复制到子代 wq；并行版将父代 wq 写入子代的 w 输入数组。这是因为并行版的递归调用 vlz_p 使用 w（第 4 个参数）作为通量输入，而串行版使用 wq。注释 correction bug le 15mai2015 标记了这是一个 bug 修正。

通量散度同步

! 串行：父子递归调用后直接计算散度

! 并行：
!$OMP BARRIER                    ! 额外同步
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
! 散度更新
!$OMP END DO NOWAIT

并行版在父子递归调用后、散度更新前插入 OMP BARRIER，确保所有线程的 wq 已就绪。注释 On rajoute ici une barriere car on veut etre sur que tous les wq soient synchronises 解释了原因。

minmaxq 对照

! 串行：
subroutine minmaxq(zq, qmin, qmax, comment)
real zq(ip1jmp1, llm)
real zzq(ip1jmp1, llm)        ! 用于打印的备份

! 并行：
subroutine minmaxq_p(zq, qmin, qmax, comment)
real zq(ip1jmp1, llm)
real zzq(iip1, jjp1, llm)     ! 三维备份，便于定位 (i,j,l)

并行版的 zzq 使用三维索引 (iip1, jjp1, llm) 而非一维展开，使得极值位置可以用 (imin, jmin, lmin) 直接定位。串行版需要从一维索引反算三维坐标。两者都受 #ifdef isminismax 保护，默认不编译。

数值等价性分析

严格等价的方面

Van Leer 斜率限制：标准斜率和乘积斜率的公式完全相同，包括 pente_max 阈值判断。
上游通量公式：q + 0.5*(1-Courant)*dq 形式一致。
Semi-Lagrangian CFL > 1 累积：逐格点追踪逻辑完全相同。
通量缩放：中间步半通量、末步全通量的策略一致。
父子 tracer 比例输送：masseq/Ratio/qperemin/masseqmin 阈值保护一致。
E-W 周期包裹：使用相同的 mod 公式。
极点 SSUM 累积：公式相同，只是加了进程条件。

可能产生差异的方面

浮点运算顺序：并行版的 OMP 循环可能以不同顺序累加，在 pente_max < -1e-5（乘积斜率）的除法中可能产生微小差异。
极点 sin/cos 投影：串行版无条件做两个极点的投影，并行版各进程只做一个。如果投影结果因浮点累加顺序不同而有微小差异，可能导致归零前的 dyq 不同——但由于最终都归零，实际无影响。
vlz 父子通量传递：w(ij,l,iq2) = wq(ij,l,iq)（并行）vs wq(ij,l,iq2) = wq(ij,l,iq)（串行）。这是代码结构差异而非数值差异，因为递归调用的入口参数不同但传入的值相同。
mw 维度扩展：mw(:,:,iq) vs mw(:,:) 是存储差异，不影响数值结果。

结论

在相同输入、1 个 MPI 进程和 1 个 OMP 线程条件下，源码结构预期 vlsplt 和 vlsplt_p 应给出相同数值结果，但本页未提供回归测试证明“位级相同”。多线程 OMP 或多进程 MPI 下，halo 交换本身传递原值而非近似值；实际是否 bitwise reproducible 仍应以串并行数值回归为准。

复现要点

维护两套代码时，任何对 Van Leer 数值格式的修改必须同时应用到 vlsplt.F 和 vlsplt_p.F。遗漏任一侧会导致串并行结果不一致。
vlsplt_p.F 中的 correction bug le 15mai2015（w vs wq）是已知的串并行不一致修正。检查类似 bug 时应关注父子 tracer 通量传递路径。
GOTO 8888 在串行版中保留但被无条件跳过，在并行版中被注释掉。两者效果相同（极点 dyq 归零），但代码维护时应以并行版的条件处理为准。
并行版的 OMP BARRIER 数量（vlz_p 中至少 4 处）是正确性关键。减少 BARRIER 前必须分析数据依赖。
mw(:,:,iq) vs mw(:,:) 的维度差异反映了串行"主例程循环 iq"和并行"每个 iq 独立调用"的调用模式差异。

待确认

Mars 并行构建是否总是使用 1 个 OMP 线程（即 OMP_CHUNK 的实际值），使得串并行数值差异在实践中不出现。
mw(ip1jmp1,llm+1,nqtot) 的内存开销在高 tracer 数配置下是否显著。
串行版 GOTO 8888 是否会在未来清理中被删除。