序

早就听说设计模式在软件领域非常重要。是时候好好学习一下了。

通常所见的设计模式教程面向的是IT软件开发者，很少有见到专为科学计算软件开发者写的教程。与企业软件相比，科学计算软件当然有着不同的特点。因此，有必要专门出一篇面向科学计算软件开发者的教程，希望能对相关从业者有帮助。我在ChatGPT的帮助下完成了这篇文章。我也是边写边学，受益颇多。

LAMMPS和SPHinXsys是两个非常优秀的大型科学计算软件，我都详细研究过源码。因此，我以这两个库为例，说明科学计算软件中的设计模式。当然，本文希望面向的并不只是熟悉这两个软件的开发者——这些设计思想可以迁移到多种科学计算领域，本文自始至终也会强调这一点，从粒子方法延伸到网格方法（有限元、有限体积）、线性求解器、优化器等其他科学计算领域。

本文字数约一万两千字，读起来当然不会容易，可以收藏起来慢慢读。每一节后都设置了一道思考题，文末给出了思考题的答案。如果您能坚持读完，我相信一定会有所收获。感兴趣的读者，甚至可以去相关软件的GitHub仓库翻阅源码，进一步了解。

1. 从科学计算中的重复结构开始

不同科学计算方法使用不同的数学语言。有限元程序遍历单元和积分点，有限体积程序遍历控制体及其界面，粒子程序遍历中心粒子和邻居粒子，稀疏矩阵算法则遍历矩阵行及其非零元素。

这些程序在数学上可能相去甚远，在软件结构上却经常具有一个共同特点：外层存在相对稳定的遍历过程，内层执行随物理模型和数值方法变化的局部计算。

以离散粒子程序为例，其基本计算过程可以简化为：

for (std::size_t i = 0; i < number_of_particles; ++i){    initialize(i);for (std::size_t j : neighbors(i))    {        interact(i, j);    }    update(i);}

neighbors(i) 给出中心粒子 i 的邻居集合，interact(i, j) 计算邻居粒子 j 对中心粒子的贡献，update(i) 则根据累积结果更新粒子的状态。

LAMMPS 和 SPHinXsys 都大量使用这种邻域计算结构。LAMMPS 通常通过邻居表组织一定截断范围内的粒子相互作用；SPHinXsys 则通过粒子关系组织核函数支撑域内的粒子求和。两者的局部公式和物理含义并不相同，但都需要回答“邻居是谁”“怎样计算贡献”“何时更新状态”这些软件设计问题。（pair_style command — LAMMPS documentation）

这里不能把 LAMMPS 简单理解为只能进行原子或分子动力学模拟。LAMMPS 官方将其描述为既能模拟原子系统，也能作为并行粒子模拟器处理介观乃至连续尺度问题。其可选 SPH package 提供传统 SPH 相关的 pair、fix 和 compute style，RHEO package 则提供更现代、灵活的混合 SPH 实现，可用于流动、多相材料和流固系统。（8.1. Package details — LAMMPS documentation）

SPHinXsys 以 SPH 为统一基础，面向流体、固体、流固耦合、扩散反应和多体动力学等多物理问题。它同样不是某一个单独流体求解器，而是一套用于建立领域应用的 C++ 科学计算库。

本文选择 LAMMPS 和 SPHinXsys，并不是要比较两种 SPH 实现，也不是把受众限制在粒子方法开发者。粒子邻域循环只是一个容易观察的软件结构。类似问题同样存在于有限元、有限体积、谱方法、离散元、线性求解器和优化程序中：

• 局部物理模型需要替换；
• 对象需要根据配置创建；
• 计算步骤具有稳定顺序；
• 多个过程需要组合；
• 高层配置需要灵活，底层内核又需要高性能。

设计模式就是用于组织这些变化的常见方法。

设计模式不是某段必须照抄的代码，也不是某一种编程语言功能。它描述的是一类反复出现的设计问题，以及对象、算法和依赖关系可以怎样安排。继承、虚函数、模板、函数对象和注册表只是实现模式的工具。

本文将依次讨论策略模式、工厂模式、模板方法模式、命令模式和组合模式。每一种模式都先从一般问题出发，再观察它在 LAMMPS 和 SPHinXsys 中的具体表现。

思考题 1

既然许多科学计算程序都可以概括为“遍历数据并执行局部计算”，为什么不能把所有算法都放入一个通用循环，再通过大量 if 或 switch 选择具体公式？

2. 策略模式：把变化的算法放到稳定接口之后

科学计算程序经常需要完成一项稳定任务，但完成任务的具体算法可能变化。

例如，程序始终需要计算粒子相互作用，但可以使用 Lennard-Jones 势、Morse 势、SPH 压力项或其他模型。有限体积程序始终需要计算界面通量，但可以选择 Roe、HLL 或 HLLC 等不同方法。非线性求解器始终需要计算搜索方向，但可以使用不同线性求解策略。

如果调用者直接包含所有算法的判断，代码可能写成：

voidcompute_interaction(    ParticleData &particles,const NeighborList &neighbors,    Model model){switch (model)    {case Model::LennardJones:        compute_lennard_jones(particles, neighbors);break;case Model::Morse:        compute_morse(particles, neighbors);break;case Model::SphPressure:        compute_sph_pressure(particles, neighbors);break;    }}

模型较少时，这种写法简单直接。模型不断增加后，调用者会逐渐依赖所有具体算法。每增加一种模型，都要修改已有分支；不同模型需要的参数、邻居关系和辅助变量也会集中到同一个函数中。

策略模式用于处理这种“任务稳定、算法可替换”的问题。

策略模式通常包含三个角色：

1. 策略接口：规定一类算法能够执行什么操作；
2. 具体策略：实现一种具体算法；
3. 上下文：使用策略完成自己的工作，但不依赖某个具体策略。

一个简化的 C++ 策略可以写成：

classInteractionStrategy{public:virtual ~InteractionStrategy() = default;virtualvoidcompute(        ParticleData &particles,const NeighborList &neighbors)= 0;};classLennardJonesInteraction    :public InteractionStrategy{public:voidcompute(        ParticleData &particles,const NeighborList &neighbors)override{// 遍历粒子和邻居，计算 Lennard-Jones 作用。    }};classSolver{public:explicitSolver(std::unique_ptr<InteractionStrategy> strategy)        : strategy_(std::move(strategy)){    }voidadvance(){        strategy_->compute(particles_, neighbors_);    }private:    ParticleData particles_;    NeighborList neighbors_;std::unique_ptr<InteractionStrategy> strategy_;};

Solver 只知道自己需要一个“粒子相互作用策略”，并不知道实际使用哪种公式。新增模型时，可以添加新的策略类，而不必修改 Solver 的时间推进流程。

这里最重要的不是虚函数，而是依赖方向。稳定的求解流程依赖抽象策略，具体物理模型实现抽象策略。

策略模式也不要求一定使用继承。如果算法在编译期已经确定，可以把策略作为模板参数；如果只有少数固定类型，可以使用 std::variant；如果算法只是一个简单公式，也可以使用函数对象。

LAMMPS 中的 pair_style

LAMMPS 将大量非键相互作用模型组织为不同的 pair_style。负责这类计算的具体类派生自 Pair 基类。新的 pair style 通常实现 compute()、settings() 和 coeff() 等接口，并可以根据需要声明邻居表、通信和重启方面的要求。（3.6. Pair styles — LAMMPS documentation）

用户可以在输入脚本中选择具体相互作用模型：

pair_style lj/cut 2.5pair_coeff * * 1.0 1.0

输入命令选择 lj/cut 后，LAMMPS 的时间推进流程不必增加“当前模型是否为 Lennard-Jones”的条件判断。核心流程只需要通过 Pair 接口调用当前模型。

从策略模式的角度看：

• Pair 是策略接口；
• 不同 pair style 是具体策略；
• LAMMPS 的力计算流程是策略使用者。

LAMMPS 的 SPH package 也沿用同一 style 体系。例如，sph/rhosum pair style 使用核函数插值计算 SPH 粒子的局部质量密度。它虽然与原子势具有不同物理含义，却仍然可以通过 LAMMPS 的 pair style 扩展机制接入计算。

这说明设计模式可以复用软件结构，却不会消除数值模型之间的区别。原子势、SPH 密度和其他粒子作用可以具有相似的调用接口，但它们对邻居表、状态变量和更新规则的要求仍然不同。

SPHinXsys 中的模板策略

SPHinXsys 也需要替换局部算法，但它的底层 C++ 接口经常使用模板参数表达具体策略。

在官方算例中，可以看到以下形式：

SimpleDynamics<NormalDirectionFromBodyShape>update_normal(body);InteractionWithUpdate<LinearGradientCorrectionMatrixInner>update_correction(inner_relation);Dynamics1Level<solid_dynamics::Integration1stHalfPK2>stress_relaxation_first_half(inner_relation);ReduceDynamics<solid_dynamics::AcousticTimeStep>estimate_time_step(body);

body：一组具有共同物理属性的粒子

dynamics：封装一次物理或数值操作的计算对象

这些声明来自 SPHinXsys 算例中反复出现的组织方式：具体物理算法作为模板参数传入通用 dynamics 类型。(SPHinXsys/tests/3d_examples/test_3d_twisting_column)

以 ReduceDynamics<solid_dynamics::AcousticTimeStep> 为例，AcousticTimeStep 规定怎样根据单个粒子的状态估计局部时间步，ReduceDynamics 则负责遍历粒子并把局部结果归约为一个全局值。

这里的具体策略在编译期确定。ReduceDynamics<AlgorithmA> 与 ReduceDynamics<AlgorithmB> 是两个不同的 C++ 类型，编译器可以针对具体组合生成代码。

因此，LAMMPS 和 SPHinXsys 都使用了策略思想，但采用了不同实现：

• LAMMPS 的很多策略通过基类和虚函数在运行时选择；
• SPHinXsys 的很多底层策略通过模板参数在编译时选择。

这不是两个不同的设计模式，而是同一模式的两种实现。

策略应该提取到哪一层

假设两个压力松弛算法只有黎曼求解器不同。可以把整个压力松弛过程设计成两个策略，也可以只把局部通量计算提取成策略。

如果把整个过程都复制到两个类中，邻居遍历和状态更新会产生重复。若只提取真正变化的局部公式，公共执行过程可以继续复用。

反过来，如果两个算法连数据依赖和更新时机都不相同，强行把它们放在同一个细粒度接口之后，会产生大量条件判断和空操作。

策略边界应当围绕已经观察到的变化确定，而不是围绕模式名称确定。

思考题 2

如果一个虚函数每个时间步只调用一次，而虚函数内部完成数百万次粒子计算，它会成为主要性能问题吗？如果虚函数在每一对邻居粒子的计算中调用，答案会发生什么变化？

3. 工厂模式：把对象的创建与使用分开

使用策略模式之后，调用者可以通过统一接口使用不同算法，但程序仍然需要决定创建哪一种具体策略。

直接创建对象时，代码可能写成：

std::unique_ptr<InteractionStrategy> strategy;if (name == "lj"){    strategy =std::make_unique<LennardJonesInteraction>();}elseif (name == "morse"){    strategy =std::make_unique<MorseInteraction>();}

这段代码把对象创建集中到了一个位置。随着具体类型不断增加，创建函数仍然会形成大型条件分支，并依赖所有具体类型。

工厂模式用于把对象创建与对象使用分开。

调用者向工厂说明需要哪类对象，工厂负责选择具体类型、执行构造，并把符合公共接口的对象返回给调用者。对象使用者不必了解具体类型的类名和构造过程。

一个简单工厂可以写成：

std::unique_ptr<InteractionStrategy>create_interaction(conststd::string &name){if (name == "lj")    {returnstd::make_unique<            LennardJonesInteraction>();    }if (name == "morse")    {returnstd::make_unique<            MorseInteraction>();    }throwstd::invalid_argument("Unknown interaction model: " + name);}

这种简单工厂已经能把创建逻辑从求解器中移出，但它本身仍然需要随着具体类型增加而修改。大型框架因而常在工厂之外再引入注册机制。

本文使用广义的“工厂”概念，重点讨论根据标识创建具体对象的机制，不进一步区分各类工厂模式。

LAMMPS 的 style 注册和工厂

LAMMPS用户通常用输入脚本写命令，而不是直接依赖某个内部 style 类。LAMMPS 根据关键字查找相应创建函数，并创建派生自公共基类的具体对象。官方开发文档把这一机制明确描述为 Factory。（4.8. Writing new styles — LAMMPS documentation）

我之前写过一篇文章LAMMPS如何根据指定style调用相应的派生类？，专门分析了这里的style注册和使用机制。

一个 pair style 可以通过类似下面的声明注册：

PairStyle(foo, PairFoo);

注册信息把用户可见的 foo 与 C++ 类型 PairFoo 联系起来。构建系统收集这些声明后，输入脚本便可以使用：

pair_style foo

这种设计避免在核心代码中手工维护包含全部 style 的条件分支。新的 style 主要负责提供自己的实现和注册信息。

LAMMPS 对 pair、fix、compute、region 和其他大量用户可选功能采用相似的扩展方式。可选 package 则把一组相关 style 和功能组织在一起，使构建时能够按需要启用。

这里可以看到三种机制的分工：

• 基类提供稳定接口；
• 具体派生类实现算法；
• 工厂和注册表根据输入创建对象。

继承本身并不能完成运行时配置。真正把用户命令与具体对象连接起来的是工厂和注册机制。

SPHinXsys 为什么较少依赖字符串工厂

SPHinXsys 的典型应用直接使用 C++ 声明具体 dynamics 类型：

Dynamics1Level<PressureRelaxation>pressure_relaxation(relation);

relation：粒子之间的邻域拓扑

类型选择已经写入程序，因此不必先把字符串转换为 C++ 类型。对于这种在编译前已经确定的算法组合，直接构造对象比建立运行时注册表更简单。

这不表示 SPHinXsys 不能使用工厂。如果应用希望从配置文件中选择几种已经编译好的算法，可以在高层增加一个工厂，再让返回的具体对象内部使用模板化内核。

补充说明：SPHinXsys 项目目前还在开发高层接口 SPHinXsim，希望通过 Python 和结构化配置降低应用构建门槛；由于它尚未正式推出，本文不把它作为主要分析对象。项目 README 已出现其 JSON 配置和验证流程的相关说明。

什么时候不需要工厂

并不是所有科学计算项目都需要运行时工厂。

如果程序只有两三种算法，并且算法在编译前已经确定，直接构造具体对象通常更加清楚。加入字符串注册表会额外引入：

• 名称管理；
• 配置错误处理；
• 对象生命周期管理；
• 注册顺序；
• 未知类型检查。

只有当用户需要通过配置选择类型、第三方需要独立添加实现，或者对象类型确实只能在运行时确定时，工厂才具有明显价值。

思考题 3

一个研究程序只有三种积分器，所有算例都在 C++ 中明确指定积分器类型。是否应该为了“可扩展性”加入字符串注册表和运行时工厂？

4. 模板方法模式：固定算法骨架，开放局部步骤

策略模式关心“某一步使用哪一种算法”。科学计算中还存在另一类问题：一个完整算法由哪些步骤组成，这些步骤应该按什么顺序执行？

许多数值算法具有相对稳定的流程。例如，一个邻域计算可能包含：

初始化局部变量  ↓遍历邻居并累积贡献  ↓更新中心粒子状态

不同模型会改变每一步的具体公式，但整体顺序可能保持不变。

模板方法模式用于处理“整体流程稳定、局部步骤变化”的问题。它由一个高层方法规定算法骨架，再把其中部分步骤交给具体实现。

一个简化的模板方法可以写成：

classParticleOperation{public:voidexecute(        ParticleData &particles,const NeighborList &neighbors){for (std::size_t i = 0;             i < particles.size(); ++i)        {            initialize(i);for (std::size_t j : neighbors(i))            {                interact(i, j);            }            update(i);        }    }protected:virtualvoidinitialize(std::size_t i)= 0;virtualvoidinteract(std::size_t i,std::size_t j)= 0;virtualvoidupdate(std::size_t i)= 0;};

execute() 规定算法顺序，派生类提供局部步骤。

需要注意，“模板方法”中的“模板”指算法流程模板，并不专指 C++ 的 template。经典模板方法常用继承和虚函数实现，也可以使用模板参数、函数对象或组合实现同样的设计意图。

LAMMPS 时间推进中的稳定骨架

LAMMPS 的一个时间步不仅包括粒子相互作用，还可能包括：

• 时间积分；
• 周期边界处理；
• 原子迁移和通信；
• 邻居表重建；
• 力计算；
• 约束和外力；
• 诊断与输出。

核心时间推进流程必须保持稳定，但用户又需要加入新的温控、约束、积分和统计方法。

LAMMPS 的 Fix 机制为时间步提供多个扩展位置。一个 fix 通过 setmask() 声明自己参与哪些阶段，并按需实现 initial_integrate()、pre_force()、post_force()、final_integrate() 或 end_of_step() 等方法。（3.9. Fix styles — LAMMPS documentation）

简化后的流程可以表示为：

initial_integrate()  ↓通信与邻居表更新  ↓pre_force()  ↓计算粒子相互作用  ↓post_force()  ↓final_integrate()  ↓end_of_step()

例如：

• 时间积分型 fix 可以参与 initial_integrate() 和 final_integrate()；
• 修改粒子受力的 fix 可以参与 post_force()；
• 统计或采样操作可以在 end_of_step() 执行。

时间步的整体顺序由 LAMMPS 控制，各个 fix 只实现自己关心的阶段。

从模板方法的角度看，LAMMPS 时间步是稳定骨架，fix 提供可变化步骤。从事件回调的角度看，fix 注册自己感兴趣的调用位置。这两个解释可以同时成立。

SPHinXsys 的 dynamics 骨架

SPHinXsys 使用多种 dynamics 包装器表达不同的粒子计算流程。官方算例中常见：

• SimpleDynamics<...>；
• InteractionWithUpdate<...>；
• Dynamics1Level<...>；
• ReduceDynamics<...>。

例如，二维溃坝算例使用：

• SimpleDynamics 施加重力和更新壁面法向；
• Dynamics1Level 进行压力和密度松弛；
• InteractionWithUpdate 进行复杂关系上的密度求和；
• ReduceDynamics 估计平流和声学时间步。

这些包装器表示的不是同一种循环。

SimpleDynamics 适合主要依赖单个粒子的局部操作；InteractionWithUpdate 表示先计算粒子相互作用，再执行状态更新；ReduceDynamics 则遍历粒子并把局部结果归约为一个整体结果。

具体物理算法作为模板参数填入这些骨架。整体执行方式稳定，局部物理公式可以变化，这正是模板方法的设计意图。

不要强迫所有算法共享一个骨架

并非所有粒子算法都严格遵循“初始化—交互—更新”。

某些算法需要在处理每个邻居后立即修改状态；某些多体势还要在邻居循环中继续遍历第三个粒子；某些隐式算法需要先装配全局系统，再交给线性求解器处理。

若为了复用一个模板方法而加入大量空函数、标志位和条件分支，说明这个所谓的公共骨架并不稳定。

正确做法可能是建立两个算法族，各自拥有自己的模板方法，而不是把所有模型压进一个万能基类。

思考题 4

某个新算法必须在处理每个邻居之后立即更新中心状态，而现有模板方法规定先完成全部邻居累积，再统一更新。应该怎样处理这个差异？

5. 命令模式：把一次计算封装成可执行对象

科学计算中的一个操作往往不只是一个公式。它可能还需要保存：

• 作用的数据对象；
• 邻居关系；
• 材料参数；
• 临时变量；
• 时间步参数；
• 并行执行策略。

如果每次调用都传递全部信息，主时间循环会充满长参数列表，算法的高层含义也会被数据细节淹没。

命令模式把一次操作及其执行所需的上下文封装为对象。调用者不必了解操作内部怎样完成，只需要触发命令。

一个简化命令接口可以写成：

classSimulationCommand{public:virtual ~SimulationCommand() = default;virtualvoidexecute(double dt)= 0;};

具体命令在构造时绑定数据：

classDensityUpdateCommand    :public SimulationCommand{public:    DensityUpdateCommand(        ParticleData &particles,        NeighborRelation &relation)        : particles_(particles),          relation_(relation)    {    }voidexecute(double dt)override{// 使用 particles_ 和 relation_ 更新密度。    }private:    ParticleData &particles_;    NeighborRelation &relation_;};

主循环只负责编排命令：

apply_gravity.execute(dt);update_density.execute(dt);compute_pressure.execute(dt);integrate_state.execute(dt);

SPHinXsys 中的 dynamics 对象

SPHinXsys 的 dynamics 对象具有明显的命令特征。对象在构造时绑定 body、relation 和模型参数，在时间循环中通过 exec() 执行。

二维溃坝算例会构造重力、密度求和、压力松弛、密度松弛和时间步估计等对象，再按照数值算法的顺序调用它们。

其调用过程可以简化为：

constant_gravity.exec();fluid_density_by_summation.exec();double dt =    fluid_acoustic_time_step.exec();fluid_pressure_relaxation.exec(dt);fluid_density_relaxation.exec(dt);

这段主循环接近数值方法的书面描述。读者可以直接看到“施加重力”“更新密度”“计算时间步”和“执行松弛”，而不必在同一层阅读所有粒子数组与邻居循环。

命令对象还使初始化与执行分离。一些 dynamics 在构造时查找变量、建立数据访问关系或准备临时状态；进入时间循环后，exec() 只完成高频计算。

LAMMPS 中的命令式对象

LAMMPS 的输入脚本本身由一系列领域命令组成，例如：

pair_style lj/cut 2.5fix integrate all nverun 10000

输入命令会创建或配置长期存在的 pair、fix 和 compute 对象，也会触发读取数据、创建原子和开始运行等一次性操作。

严格来说，LAMMPS 的所有对象都不能简单归入经典命令模式。例如，Pair 更适合作为策略，Fix 同时具有回调和插件特征。但从高层看，LAMMPS 输入系统确实把用户意图封装为一系列可解析和执行的领域操作。

这提醒我们：成熟框架中的一个对象可能同时体现多种模式。判断模式时应关注它在当前问题中承担的职责，而不是试图给每个类贴上唯一标签。

命令模式与策略模式的区别

策略和命令都可以把行为封装成对象，但关注点不同。

策略模式回答“完成同一任务可以采用哪一种算法”。例如，为压力项选择不同离散公式。

命令模式回答“需要执行哪一项操作”。例如，更新密度、施加边界条件或输出结果。

策略通常被某个上下文长期持有，用于改变上下文的算法；命令则更强调一个可安排、可执行的动作。

一个命令内部完全可以使用某个策略。例如，“更新压力”命令可以保存一个黎曼求解器策略。

命令对象不能隐藏数据依赖

把所有操作都封装为 exec() 后，主循环会变得整洁，但也可能掩盖步骤之间的数据关系。

例如：

• 密度更新需要读取最新邻居关系；
• 压力计算需要读取已经更新的密度；
• 粒子排序后，部分数据访问关系需要重新建立；
• 输出命令可能要求所有并行数据已经同步。

SPHinXsys 官方算例会特别提示，一些数值方法的构造顺序存在数据依赖；粒子排序、cell linked list 更新和 configuration 更新也必须保持正确顺序。

configuration：某一时刻实际建立的邻居配置

因此，统一的 exec() 接口并不表示所有命令可以任意交换。可执行性和可交换性是两件不同的事。

思考题 5

两个计算对象都提供无参数的 exec()，是否说明它们可以任意交换执行顺序？程序应该怎样表达隐藏的数据依赖？

6. 组合模式：在不同层次组织复杂计算

实际科学计算很少只包含一个物理过程。

一个粒子系统可能同时受到短程作用、长程作用、外力和边界约束；一个流固耦合问题可能同时包含流体内部作用、固体内部作用、流体与固体之间的接触，以及多个时间尺度上的积分过程。

最直接的实现方式，是为每一种物理组合编写一个新的求解器类。例如，程序可能先有 FluidSolver 和 SolidSolver，随后又出现 FluidWallSolver、FluidSolidSolver、FluidSolidWallSolver。随着可选过程增加，组合数量会快速增长，大量已有算法也会在新类中被重复实现。

组合模式用于解决这类问题。

组合模式的基本思想是：用较小的对象构造更复杂的对象，并尽量让调用者以相似的方式使用单个对象和组合对象。

经典组合模式通常包含两类对象：

1. 叶子对象：完成一项基本操作；
2. 组合对象：内部保存多个叶子对象或其他组合对象，并把操作转发给它们。

一个简化实现可以写成：

classInteraction{public:virtual ~Interaction() = default;virtualvoidcompute()= 0;};classCompositeInteraction :public Interaction{public:voidadd(std::unique_ptr<Interaction> interaction){        interactions_.push_back(std::move(interaction));    }voidcompute()override{for (auto &interaction : interactions_)        {            interaction->compute();        }    }private:std::vector<std::unique_ptr<Interaction>>        interactions_;};

单个相互作用和多个相互作用构成的组合都提供 compute()。调用者不必知道自己操作的是一个基本对象，还是包含多个子过程的组合对象。

LAMMPS 的 pair_style hybrid

LAMMPS 也提供了明确的行为组合机制。pair_style hybrid、hybrid/overlay 和 hybrid/scaled 可以在同一个模拟中组织多个 pair style。

hybrid 通常让不同原子类型对使用不同的子 style；hybrid/overlay 允许同一类型对同时采用多个子 style，并叠加它们的贡献；hybrid/scaled 则进一步为各个子 style 设置缩放因子。

例如，可以组合短程 Lennard-Jones 作用和库仑作用：

pair_style hybrid/overlay lj/cut 2.5 coul/long 2.0pair_coeff * * lj/cut 1.0 1.0pair_coeff * * coul/long

调用者面对的是一个 hybrid pair style，而 hybrid style 内部管理多个具体 pair style。它把若干相互作用策略组合成一个更复杂的策略，因此也具有组合模式的特征。

LAMMPS 还允许一个时间步中同时存在多个 fix。一个 fix 可以负责积分，一个 fix 施加约束，另一个 fix 完成温控或统计。它们共同构成完整模拟行为。

不过，多个 fix 按时间步阶段执行，更接近插件或工作流组合；pair_style hybrid 则更接近“组合对象内部包含多个同类子策略”的经典组合结构。

SPHinXsys的ComplexRelation和ComplexInteraction

科学计算中的组合并不只发生在物理模型层。邻域拓扑、局部相互作用和完整时间推进流程都可能需要组合。SPHinXsys 的 ComplexRelation 和 ComplexInteraction 正好展示了两个不同层次的组合。

ComplexRelation：组合邻域拓扑

在 SPHinXsys 中，一个 body 的计算既可能依赖自身内部粒子，也可能依赖其他 body 的粒子。

例如，流体粒子的压力计算可能同时需要：

• 同一流体 body 中的邻居粒子；
• 壁面 body 中的邻居粒子；
• 柔性结构 body 中的邻居粒子。

这些关系分别具有不同含义。内部粒子关系通常由 InnerRelation 表示，不同 body 之间的关系通常由 ContactRelation 表示。

如果主程序必须分别管理并更新所有关系，代码可能写成：

fluid_inner_relation.updateConfiguration();fluid_wall_contact.updateConfiguration();fluid_solid_contact.updateConfiguration();

随着接触对象增加，调用者需要了解越来越多的关系对象。

ComplexRelation 将一个内部关系和一个或多个接触关系保存到同一个对象中。源码中，它继承自 SPHRelation，内部保存一个 BaseInnerRelation 引用和一组 BaseContactRelation 指针，并提供统一的 updateConfiguration()。

它的结构可以简化为：

classComplexRelation :public SPHRelation{public:    ComplexRelation(        BaseInnerRelation &inner_relation,std::vector<BaseContactRelation *>            contact_relations);voidupdateConfiguration()override{        inner_relation_.updateConfiguration();for (size_t k = 0; k != contact_relations_.size(); ++k)            contact_relations_[k]->updateConfiguration();    }private:    BaseInnerRelation &inner_relation_;std::vector<BaseContactRelation *>        contact_relations_;};

从组合模式的角度看：

• inner relation 和 contact relation 是基本关系对象；
• ComplexRelation 是组合关系对象；
• 组合对象统一管理多个邻域 configuration；
• 调用者可以通过一个接口更新整个关系组合。

因此，ComplexRelation 可以理解为结构层面的组合。它组合的不是物理公式，而是“谁与谁可能发生作用”的邻域拓扑。

不过，不能把它解释得过于宽泛。源码注释说明，当前 ComplexRelation 主要用于同时更新多个 configuration；构造局部 dynamics 时，通常仍然分别使用 inner relation 和 contact relation。

更准确的说法是：

ComplexRelation 将内部关系和多个接触关系组合起来，使这些关系能够作为一个整体更新，但具体局部算法仍可以分别读取相应的基本关系。

这种设计也说明，组合对象不一定完全取代基本对象。它可以只在某个管理层次提供统一操作，而把更细粒度的对象继续暴露给具体算法。

ComplexInteraction：组合局部相互作用

ComplexRelation 组合的是邻域拓扑，ComplexInteraction 组合的则是作用于中心粒子上的局部算法。

一个中心粒子可能同时需要计算多类邻居贡献。例如，流体粒子可能同时受到：

• 流体内部粒子的压力与黏性作用；
• 壁面粒子的边界作用；
• 固体粒子的接触或耦合作用。

一种写法是在主循环中分别调用这些 interaction：

inner_interaction.interaction(i, dt);wall_interaction.interaction(i, dt);solid_interaction.interaction(i, dt);

这种写法可以工作，但外层执行骨架必须知道一个具体物理过程由多少种 interaction 组成。

ComplexInteraction 将多个局部 interaction 组合成一个整体。源码将其描述为“集成多个局部 dynamics”的类，典型内容包括一个 inner interaction、一个或多个 contact interaction，以及边界条件。

它没有使用运行时的对象容器，而是通过可变参数模板和递归组合多个类型。简化后的结构可以表示为：

template<classFirst, class... Others>classComplexInteraction    :public First{public:voidinteraction(std::size_t i,double dt){        First::interaction(i, dt);        others_.interaction(i, dt);    }private:    ComplexInteraction<Others...> others_;};

实际源码中的 interaction() 先执行第一个局部 interaction，再递归调用其余 interaction。对外部执行骨架来说，这组局部算法表现为一个统一的 interaction()。

从组合模式的角度看：

• 单个 inner interaction 或 contact interaction 是基本算法；
• ComplexInteraction 是组合算法；
• 组合对象内部包含其他 interaction；
• 外层循环只调用一次统一的 interaction()。

这比在主循环中手工排列多个 interaction 更接近经典组合模式。

与经典组合模式写法不同的是，ComplexInteraction 的组合关系主要在编译期建立。子对象的类型由模板参数决定，组合结构通过递归模板生成，而不是在运行时向 std::vector 中加入基类指针。

因此，从设计意图上，可以把 ComplexInteraction 称为组合模式的静态实现。它体现了组合模式的设计意图，却使用 C++ 模板而不是虚函数和运行时对象树。

这再次说明，设计模式不等于某一种固定语法。组合模式既可以通过运行时多态实现，也可以通过模板和编译期类型组合实现。

工作流组合不等于经典组合模式

在科学计算程序中，主循环经常写成：

apply_gravity.exec();update_density.exec();compute_viscosity.exec();pressure_relaxation.exec(dt);density_relaxation.exec(dt);

这些命令共同构成完整算法，因此也可以说它们形成了一种行为组合。

但是，它们通常没有共同的组合对象，也不一定共享完全相同的接口和返回值。调用者仍然明确知道每一个步骤，并负责维护执行顺序。

因此，更准确的分类是：

• pair_style hybrid：相互作用策略的运行时组合；
• ComplexRelation：结构组合；
• ComplexInteraction：局部算法的静态组合；
• 多个 exec() 或 fix 的排列：工作流组合。

软件可组合不等于物理可组合

组合模式解决的是软件结构问题：怎样把多个小对象组织成较复杂的对象。

它不能回答这些对象在物理和数值上能否直接组合。

在组合两个相互作用模型前，仍然需要检查：

• 是否重复计算了同一种物理作用；
• 两个贡献是否具有相同单位；
• 能量、力或通量是否可以直接相加；
• 是否满足质量、动量或能量守恒；
• 两个算法是否需要不同邻居范围；
• 时间步稳定条件是否兼容；
• 更新顺序是否影响数值结果；
• 多体项是否在拆分后仍然保持完整。

ComplexInteraction 也存在类似问题。模板能够把多个 interaction 编译到同一个局部计算中，却不能自动判断它们的执行顺序是否合理，也不能证明多个贡献满足所需守恒关系。

因此，科学计算中的组合应当分成两个判断：

1. 软件接口是否允许组合；
2. 数学和物理模型是否允许组合。

只有两个问题的答案都为肯定，组合才真正成立。

思考题 6

ComplexRelation、ComplexInteraction 和多个 dynamics 的顺序调用都体现了组合思想，但三者组合的对象和提供的统一接口并不相同。为什么区分结构组合、局部算法组合和工作流组合很重要？如果把它们全部放入一个通用的 Composite 类，会产生哪些问题？

7. 运行时多态、静态多态与科学计算性能

LAMMPS 和 SPHinXsys 经常被概括为两种不同风格：

• LAMMPS 主要使用继承和运行时多态；
• SPHinXsys 主要使用模板和静态多态。

这个概括能够帮助入门，但不应被理解为非此即彼。

设计模式描述的是算法和对象之间的关系。动态多态和静态多态只是实现关系的不同方式。

运行时多态

运行时多态通常使用基类指针和虚函数。实际调用哪一个实现，由对象在运行时的具体类型决定。

它适合以下场景：

• 用户通过输入文件选择模型；
• 类型在运行前不能确定；
• 第三方模块需要加入新的具体实现；
• 高层需要保存一组接口相同、类型不同的对象。

LAMMPS 的 style 系统符合这些需求。用户提供 style 名称，工厂创建具体对象，主流程通过公共基类调用它。

运行时多态的性能影响取决于分派位置。

下面的调用每个时间步只发生一次：

pair_model->compute();

compute() 内部可能进行数百万次粒子相互作用。与内存访问、邻居遍历和浮点计算相比，一次虚函数调用通常不是主要成本。

下面的写法则把虚函数放入最内层：

for (std::size_t i = 0;     i < particle_count; ++i){for (std::size_t j : neighbors(i))    {        pair_model->interact(i, j);    }}

此时动态分派可能发生数百万乃至更多次，也会妨碍具体局部公式的内联。问题不在于“虚函数一定慢”，而在于虚函数位于哪一个计算层级。

静态多态

静态多态通常通过模板实现。具体类型在编译时确定，编译器可以针对每一种组合生成代码。

它适合以下场景：

• 算法组合在编译前已经确定；
• 局部函数位于高频热点；
• 需要内联和常量传播；
• 需要为不同执行后端生成专用 kernel。

SPHinXsys 的 dynamics 组合体现了这种思路。具体算法作为模板参数传入执行骨架，编译器能够看到局部操作的具体类型。

静态多态也有成本：

• 类型表达更加复杂；
• 编译时间可能增加；
• 不同组合可能导致代码膨胀；
• 错误信息较难阅读；
• 未编译的新类型不能只靠配置文件接入。

因此，模板不是自动获得高性能的保证。数据布局、访存模式、循环结构和并行冲突往往比函数调用方式更重要。

粗粒度动态，细粒度静态

许多科学计算软件适合采用混合结构：

用户输入或配置  ↓运行时工厂选择高层模型  ↓创建粗粒度求解对象  ↓调用模板化计算内核  ↓遍历粒子、网格或矩阵数据

运行时多态负责高层选择，静态多态负责热点内核。用户获得配置灵活性，编译器也能够优化内层循环。

例如，可以预先编译三种核函数对应的 SPH 求解对象，再由高层工厂根据输入选择其中一个。动态分派只发生在调用整个密度计算时，具体求解对象内部的每对粒子核函数仍然可以被内联。

数据导向设计同样重要

科学计算中的设计模式不能只讨论类关系。

LAMMPS 的高层结构使用 C++ 类，而实际计算方法大量操作简单向量和数组。SPHinXsys 则采用结构化数组及粒子排序等数据组织方式，以改善缓存访问。（3.1. Overview — LAMMPS documentation）

高层对象可以表达：

• 模型；
• 材料；
• 邻居关系；
• 时间推进步骤；
• 输出与测试。

底层数据仍应尽量满足：

• 连续存储；
• 较少指针跳转；
• 良好缓存局部性；
• 方便向量化；
• 方便 GPU 合并访存；
• 清楚的并行写入规则。

面向对象设计不意味着每个粒子都应成为带虚函数和独立内存分配的复杂对象。设计模式主要管理高层变化，数据导向设计主要保证底层效率，两者应在不同层次协作。

这些原则不限于粒子方法

在有限元程序中，可以把材料本构看作策略，把单元创建看作工厂，把“形函数计算—积分点循环—装配”看作模板方法。

在有限体积程序中，可以把数值通量和限制器看作策略，把边界条件封装为命令或回调，把多种物理通量组合为复合模型。

在线性代数库中，可以把预条件器看作策略，把求解器配置转换为对象的过程看作工厂，把迭代求解的“初始化—残差—更新—收敛判断”看作模板方法。

在优化软件中，可以把目标函数、梯度方法、线搜索和停止准则分别组织为策略，再通过高层求解流程组合。

模式的名称没有变化，具体数据结构和数值含义却会随领域改变。

思考题 7

一个重构把五种算法统一到同一个模板框架中，删除了大量重复循环，但每一种算法都需要增加条件分支和空操作才能适应公共骨架。这个设计可能出现了什么问题？

8. 本课小结

LAMMPS 和 SPHinXsys 都包含大量围绕数据遍历、邻域关系和局部物理算法展开的计算。二者也都不仅限于单一传统应用：LAMMPS 可以通过 SPH 和 RHEO 等 package 扩展到 SPH 与多尺度粒子问题，SPHinXsys 则用统一框架组织流体、固体和多物理耦合。

策略模式隔离可以替换的算法，工厂模式把对象创建与使用分开，模板方法固定计算骨架，命令模式封装计算步骤，组合模式则用较小对象构造复杂系统。

这些模式并不要求使用某一种 C++ 语法。LAMMPS 较多在高层使用运行时多态和注册工厂，SPHinXsys 较多在底层使用模板和静态组合；两者同时也都依赖规则数据结构和高性能内核。

对科学计算开发人员而言，真正重要的不是记住模式类图，而是识别变化发生在哪个层次，并让高层保持可扩展、底层保持可验证和高效。

9. 思考题参考答案

问题 1

通用循环只能表达表面的遍历相似性，不能表达不同算法的数据需求、邻居形式、更新规则和并行约束。

随着模型增加，一个包含全部算法的 switch 会依赖所有具体实现。新增模型时必须修改稳定代码，模型参数和特殊情况也会逐渐集中到同一个函数中。

更重要的是，有些差异并不只发生在局部公式。某些算法使用 half neighbor list，某些算法需要 full list；某些算法只更新中心对象，某些算法还要更新邻居；某些算法需要额外遍历第三个对象。

如果只有少量简单且稳定的算法，条件分支仍可能是合理选择。只有当变化开始反复出现时，才需要引入模式。

问题 2

如果虚函数每个时间步只调用一次，而函数内部完成大量计算，一次动态分派通常不会成为主要性能问题。

如果虚函数位于每一对邻居的计算中，调用频率会大幅增加，编译器也难以内联具体公式。此时可以把动态多态提升到粗粒度层级：虚函数负责调用完整计算过程，具体策略内部使用普通函数或模板实现热点循环。

判断性能时应测量实际程序，而不是仅凭“虚函数慢”或“模板快”的印象。

问题 3

通常没有必要。

如果算法类型在编译前已经确定，直接构造对象或使用模板参数更加简单。字符串注册表会增加名称解析、错误处理和生命周期管理，却没有解决实际的运行时需求。

当用户需要通过输入文件选择算法，或第三方需要独立加入新类型时，工厂和注册机制才具有明显价值。

问题 4

不应强迫新算法适应错误的执行骨架。

“逐邻居立即更新”与“完成全部邻居累积后统一更新”具有不同的数据依赖和数值语义。可以为新算法建立另一种模板方法，同时复用仍然相同的邻居访问和数据准备部分。

公共骨架应来自真正稳定的流程，而不是为了消除代码重复而人为制造的一致性。

问题 5

不能。

统一的 exec() 只说明两个对象都可以被执行，不说明它们之间没有数据依赖。

应检查一个对象是否读取另一个对象写入的数据，是否共享临时状态，是否要求最新邻居关系，以及顺序变化是否改变时间离散。

依赖可以通过高层流程、明确的输入输出类型、任务图、状态标记或文档约定表达。对于关键依赖，仅依赖调用者记忆通常不够可靠。

问题 6

这道题的关键，是区分三种“组合”发生在哪一个层次。

ComplexRelation 组合的是邻域拓扑，回答“中心粒子需要和哪些粒子建立关系”；ComplexInteraction 组合的是局部相互作用，回答“针对这些邻居要计算哪些贡献”；多个 dynamics 的顺序调用组合的是完整计算流程，回答“一个时间步中各个数值步骤按什么顺序执行”。

区分这三层很重要，因为它们具有不同的输入、输出和约束。关系组合主要涉及邻居表与 configuration 的更新；局部算法组合主要涉及粒子贡献的累积、执行顺序和守恒性；工作流组合则涉及时间推进、数据依赖和不同时间尺度。若把三者混在一起，程序会难以判断一次变化究竟属于拓扑、物理公式还是流程调度。

如果把它们全部塞进一个通用 Composite 类，通常会出现以下问题：

• 接口会过于宽泛。统一接口可能同时要求 updateConfiguration()、interaction() 和 exec()，导致大量对象只能实现其中一部分，其余方法成为空操作。
• 类型语义会被削弱。一个 relation 和一个 dynamics 虽然都“可组合”，但前者描述数据关系，后者执行计算，它们并不是同一种职责。
• 数据依赖会被隐藏。多个 dynamics 的执行顺序可能影响结果，而 relation 的组合通常不表达时间顺序；通用组合类难以区分这两类约束。
• 编译期与运行时机制会被混淆。ComplexInteraction 适合用模板静态组合，关系对象可能使用引用和容器，工作流又可能需要运行时调度。强行统一会牺牲其中某些层次的表达能力或性能。
• 错误组合更容易发生。通用接口可能允许把一个 relation 当作计算步骤加入工作流，或者把需要不同更新阶段的 interaction 任意排列，导致接口合法但数值意义错误。

因此，更好的做法是让每个层次拥有自己的组合抽象，同时在更高层明确连接它们：

Relation：确定邻域拓扑  ↓Interaction：计算邻居贡献  ↓Dynamics：组织完整操作  ↓Time loop：安排数值流程

组合模式的目的不是把所有对象统一成同一种类型，而是在同一职责层次内把基本对象组织成复杂对象。对科学计算程序而言，保持拓扑、局部算法和时间流程之间的边界，往往比追求一个万能的组合接口更重要。

问题 7

公共抽象可能选择了错误层级。

五种算法虽然都有循环，但其初始化、相互作用或更新过程可能并不真正相同。大量条件和空操作说明公共骨架并不稳定。

可以缩小公共部分，只复用真正一致的邻居访问；也可以把算法分为多个家族，每个家族使用自己的模板方法。对于结构特殊或性能关键的算法，保留专用实现也是合理选择。

消除重复代码不是抽象的唯一目标。好的抽象还应让数值语义更清楚，而不是把差异隐藏在条件分支中。