使用开源工具的科学计算指南

从历史上看，计算科学在很大程度上仅限于研究科学家和博士生的领域。 然而，这些年来——也许更大的软件社区不知道——我们科学计算的蛋头们一直在悄悄地编译协作开源库来处理绝大多数繁重的工作。 结果是，现在实现数学模型比以往任何时候都更容易，虽然该领域可能还没有为大规模消费做好准备，但成功实施的门槛已经大大降低。 从头开始开发新的计算代码库是一项艰巨的任务，通常以数年为衡量标准，但这些开源科学计算项目使得运行可访问的示例成为可能，以相对快速地利用这些计算能力。

由于科学计算的目的是为自然界中存在的真实系统提供科学见解，因此该学科代表了使计算机软件接近现实的前沿。 为了使软件能够以非常高的精度和分辨率模拟现实世界，必须调用复杂的微分数学，这需要在大学、国家实验室或企业研发部门之外很少发现的知识。 最重要的是，当尝试使用零和一的离散语言来描述现实世界的连续、无穷小的结构时，会出现重大的数值挑战。 为了呈现计算上易于处理且呈现有意义结果的算法，需要进行仔细的数值转换的详尽努力。 换句话说，科学计算是重任。

科学计算的开源工具

我个人特别喜欢 FEniCS 项目，将它用于我的论文工作，并将证明我在选择它作为本教程的代码示例时的偏见。 （还有其他非常高质量的项目，例如 DUNE，也可以使用。）

FEniCS 自称是“一个为自动化科学计算开发创新概念和工具的合作项目，特别关注通过有限元方法自动求解微分方程。” 这是一个强大的库，用于解决大量问题和科学计算应用程序。 它的贡献者包括 Simula 研究实验室、剑桥大学、芝加哥大学、贝勒大学和 KTH 皇家理工学院，在过去十年中，他们共同将其打造为宝贵的资源（参见 FEniCS 代码群）。

令人惊讶的是 FEniCS 库为我们提供了多少努力。 要了解该项目所涵盖主题的惊人深度和广度，可以查看他们的开源教科书，其中第 21 章甚至比较了解决不可压缩流的各种有限元方案。

该项目在幕后为我们集成了大量的开源科学计算库，这些库可能会感兴趣或直接使用。 这些包括，没有特别的顺序，FEniCS项目调用的项目：

• PETSc：一套数据结构和例程，用于由偏微分方程建模的科学应用的可扩展（并行）解决方案。
• Trilinos 项目：一套强大的算法和技术，用于求解线性和非线性方程，由桑迪亚国家实验室开发。
• uBLAS：“一个 C++ 模板类库，为密集、压缩和稀疏矩阵提供 BLAS 1、2、3 级功能，以及许多用于线性代数的数值算法。”
• GMP：任意精度算术的免费库，对有符号整数、有理数和浮点数进行操作。
• UMFPACK：一组用于求解非对称稀疏线性系统的例程，Ax=b，使用非对称 MultiFrontal 方法。
• ParMETIS：一个基于 MPI 的并行库，它实现了用于划分非结构化图、网格和计算稀疏矩阵的填充减少排序的各种算法。
• NumPy：使用 Python 进行科学计算的基础包。
• CGAL：C++ 库形式的高效可靠的几何算法。
• SCOTCH：用于顺序和并行图形分区、静态映射和聚类、顺序网格和超图分区以及顺序和并行稀疏矩阵块排序的软件包和库。
• MPI：由学术界和工业界的一组研究人员设计的标准化和便携式消息传递系统，可在各种并行计算机上运行。
• VTK：用于 3D 计算机图形、图像处理和可视化的开源、免费提供的软件系统。
• SLEPc：用于解决并行计算机上的大规模稀疏特征值问题的软件库。

这个集成到项目中的外部包列表让我们了解它的继承能力。 例如，对 MPI 的集成支持允许在计算集群环境中跨远程工作人员进行扩展（即，此代码将在超级计算机或您的笔记本电脑上运行）。

值得注意的是，除了科学计算之外，还有许多应用程序可以利用这些项目，包括金融建模、图像处理、优化问题，甚至可能是视频游戏。 例如，可以创建一个视频游戏，该游戏使用其中一些开源算法和技术来解决二维流体流动问题，例如玩家将与之交互的洋流/河流（也许尝试和用不同的风和水流航行船）。

示例应用程序：利用开源进行科学计算

在这里，我将尝试通过展示如何在这些开源库之一（在本例中为 FEniCS 项目）中开发和实现基本的计算流体动力学方案来了解开发数值模型所涉及的内容。 FEnICS 提供 Python 和 C++ 的 API。 对于这个例子，我们将使用他们的 Python API。

我们将讨论一些相当技术性的内容，但目标只是让我们了解开发这种科学计算代码需要什么，以及今天的开源工具为我们做了多少繁琐的工作。 在此过程中，希望我们能帮助揭开复杂的科学计算世界的神秘面纱。 （请注意，提供了一个附录，其中详细介绍了对这些详细程度感兴趣的人的所有数学和科学基础。）

免责声明：对于那些在科学计算软件和应用程序方面很少或没有背景的读者，本示例的部分内容可能会让您感觉如下：

如果是这样，请不要绝望。 这里的主要内容是现有开源项目可以在多大程度上简化许多这些任务。

考虑到这一点，让我们从查看不可压缩 Navier-Stokes 的 FEnICS 演示开始。 此演示模拟流经 L 形弯管（如管道）的不可压缩流体的压力和速度。

链接的演示页面上的描述对运行代码的必要步骤进行了出色、简洁的设置，我鼓励您快速查看其中涉及的内容。 总而言之，该演示将求解不可压缩流动方程通过弯曲处的速度和压力。 该演示对随时间流动的流体进行了简短的模拟，并对结果进行了动画处理。 这是通过设置表示管道空间的网格，并使用有限元方法对网格上每个点的速度和压力进行数值求解来实现的。 然后我们遍历时间，更新速度和压力场，再次使用我们可以使用的方程。

该演示按原样运行良好，但我们将对其稍作修改。 该演示使用 Chorin 拆分，但我们将使用受 Kim 和 Moin 启发的略有不同的方法，我们希望这种方法更稳定。 这只需要我们改变用于近似对流项和粘性项的方程，但要做到这一点，我们需要存储前一个时间步的速度场，并在更新方程中添加两个附加项，这将使用之前的信息以获得更准确的数值近似。

所以让我们做出这个改变。 首先，我们在设置中添加一个新的Function对象。 这是一个表示抽象数学函数的对象，例如向量或标量场。 我们称它为un1它将存储先前的速度场在我们的函数空间V上。

 ... # Create functions (three distinct vector fields and a scalar field) un1 = Function(V) # the previous time step's velocity field we are adding u0 = Function(V) # the current velocity field u1 = Function(V) # the next velocity field (what's being solved for) p1 = Function(Q) # the next pressure field (what's being solved for) ...

接下来，我们需要在模拟的每个步骤中更改“暂定速度”的更新方式。 该字段表示忽略压力时下一个时间步长的近似速度（此时压力尚不清楚）。 这是我们用更新的 Kim 和 Moin 分数步法替换 Chorin 分裂法的地方。 换句话说，我们将更改字段F1的表达式：

代替：

 # Tentative velocity field (a first prediction of what the next velocity field is) # for the Chorin style split # F1 = change in the velocity field + # convective term + # diffusive term - # body force term F1 = (1/k)*inner(u - u0, v)*dx + \ inner(grad(u0)*u0, v)*dx + \ nu*inner(grad(u), grad(v))*dx - \ inner(f, v)*dx

和：

 # Tentative velocity field (a first prediction of what the next velocity field is) # for the Kim and Moin style split # F1 = change in the velocity field + # convective term + # diffusive term - # body force term F1 = (1/k)*inner(u - u0, v)*dx + \ (3.0/2.0) * inner(grad(u0)*u0, v)*dx - (1.0/2.0) * inner(grad(un1)*un1, v)*dx + \ (nu/2.0)*inner(grad(u+u0), grad(v))*dx - \ inner(f, v)*dx

以便演示现在使用我们更新的方法在使用F1时求解中间速度场。

最后，确保我们在每个迭代步骤结束时更新先前的速度场un1

 ... # Move to next time step un1.assign(u0) # copy the current velocity field into the previous velocity field u0.assign(u1) # copy the next velocity field into the current velocity field ...

因此，以下是我们的 FEniCS CFD 演示的完整代码，其中包括我们的更改：

 """This demo program solves the incompressible Navier-Stokes equations on an L-shaped domain using Kim and Moin's fractional step method.""" # Begin demo from dolfin import * # Print log messages only from the root process in parallel parameters["std_out_all_processes"] = False; # Load mesh from file mesh = Mesh("lshape.xml.gz") # Define function spaces (P2-P1) V = VectorFunctionSpace(mesh, "Lagrange", 2) Q = FunctionSpace(mesh, "Lagrange", 1) # Define trial and test functions u = TrialFunction(V) p = TrialFunction(Q) v = TestFunction(V) q = TestFunction(Q) # Set parameter values dt = 0.01 T = 3 nu = 0.01 # Define time-dependent pressure boundary condition p_in = Expression("sin(3.0*t)", t=0.0) # Define boundary conditions noslip = DirichletBC(V, (0, 0), "on_boundary && \ (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \ (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))") inflow = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS") outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS") bcu = [noslip] bcp = [inflow, outflow] # Create functions un1 = Function(V) u0 = Function(V) u1 = Function(V) p1 = Function(Q) # Define coefficients k = Constant(dt) f = Constant((0, 0)) # Tentative velocity field (a first prediction of what the next velocity field is) # for the Kim and Moin style split # F1 = change in the velocity field + # convective term + # diffusive term - # body force term F1 = (1/k)*inner(u - u0, v)*dx + \ (3.0/2.0) * inner(grad(u0)*u0, v)*dx - (1.0/2.0) * inner(grad(un1)*un1, v)*dx + \ (nu/2.0)*inner(grad(u+u0), grad(v))*dx - \ inner(f, v)*dx a1 = lhs(F1) L1 = rhs(F1) # Pressure update a2 = inner(grad(p), grad(q))*dx L2 = -(1/k)*div(u1)*q*dx # Velocity update a3 = inner(u, v)*dx L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx # Assemble matrices A1 = assemble(a1) A2 = assemble(a2) A3 = assemble(a3) # Use amg preconditioner if available prec = "amg" if has_krylov_solver_preconditioner("amg") else "default" # Create files for storing solution ufile = File("results/velocity.pvd") pfile = File("results/pressure.pvd") # Time-stepping t = dt while t < T + DOLFIN_EPS: # Update pressure boundary condition p_in.t = t # Compute tentative velocity step begin("Computing tentative velocity") b1 = assemble(L1) [bc.apply(A1, b1) for bc in bcu] solve(A1, u1.vector(), b1, "gmres", "default") end() # Pressure correction begin("Computing pressure correction") b2 = assemble(L2) [bc.apply(A2, b2) for bc in bcp] solve(A2, p1.vector(), b2, "cg", prec) end() # Velocity correction begin("Computing velocity correction") b3 = assemble(L3) [bc.apply(A3, b3) for bc in bcu] solve(A3, u1.vector(), b3, "gmres", "default") end() # Plot solution plot(p1, title="Pressure", rescale=True) plot(u1, title="Velocity", rescale=True) # Save to file ufile << u1 pfile << p1 # Move to next time step un1.assign(u0) u0.assign(u1) t += dt print "t =", t # Hold plot interactive()

运行程序会显示肘部周围的流动。 自己运行科学计算代码以查看动画！ 最后一帧的屏幕如下所示。

模拟结束时弯曲处的相对压力，按大小缩放和着色（无量纲值）：

模拟结束时弯曲处的相对速度作为矢量字形按大小（无量纲值）缩放和着色。

因此，我们所做的是采用现有的演示，它恰好很容易地实现了一个与我们非常相似的方案，并通过使用来自上一个时间步的信息对其进行了修改以使用更好的近似值。

在这一点上，您可能会认为这是一个微不足道的编辑。 确实如此，这在很大程度上是重点。 这个开源科学计算项目使我们能够通过更改四行代码快速实现修改后的数值模型。 在大型研究代码中，此类更改可能需要数月时间。

该项目还有许多其他演示可以用作起点。 甚至在该项目上构建了许多开源应用程序，它们实现了各种模型。

结论

科学计算及其应用确实很复杂。 没有办法解决这个问题。 但正如在许多领域中越来越真实的那样，可用开源工具和项目的不断增长的格局可以显着简化原本极其复杂和乏味的编程任务。 或许，科学计算变得足够容易，可以发现自己在研究界之外很容易被利用的时候已经很近了。

附录：科学和数学基础

对于那些感兴趣的人，这里是我们上面的计算流体动力学指南的技术基础。 下面的内容将作为一个非常有用和简明的主题摘要，这些主题通常涵盖在十几个研究生水平的课程中。 对深入理解该主题感兴趣的研究生和数学类型可能会发现这些材料非常吸引人。

流体力学

一般来说，“建模”是用级数逼近求解某个实际系统的过程。 该模型通常会涉及不适合计算机实现的连续方程，因此必须用数值​​方法进一步近似。

对于流体力学，让我们从基本方程 Navier-Stokes 方程开始本指南，并使用它来开发 CFD 方案。

Navier-Stokes 方程是一系列偏微分方程 (PDE)，它们很好地描述了流体流动，因此是我们的起点。 它们可以从通过雷诺传输定理并应用高斯定理并调用斯托克假设而得到的质量、动量和能量守恒定律推导出来。 这些方程需要一个连续假设，假设我们有足够的流体粒子来给出统计特性，例如温度、密度和速度含义。 此外，还需要表面应力张量和应变率张量之间的线性关系、应力张量的对称性和各向同性流体假设。 了解我们在此开发过程中所做和继承的假设非常重要，这样我们才能评估结果代码的适用性。 爱因斯坦表示法的 Navier-Stokes 方程，不用多说：

质量守恒：

动量守恒：

能量守恒：

其中偏应力为：

虽然非常普遍，可以控制物理世界中的大多数流体流动，但它们直接并没有多大用处。 这些方程的已知精确解相对较少，而 1,000,000 美元的千禧年奖颁发给任何能够解决存在性和平滑性问题的人。 重要的部分是我们通过一系列假设来降低复杂性（它们是经典物理学中一些最难的方程），从而为模型的开发提供了一个起点。

为了使事情保持“简单”，我们将使用我们的领域特定知识对流体做出不可压缩的假设，并假设恒定温度，这样就不需要（解耦）成为热方程的能量守恒方程。 我们现在有两个方程，仍然是 PDE，但在解决大量实际流体问题的同时要简单得多。

连续性方程

动量方程

在这一点上，我们现在有一个很好的不可压缩流体流动的数学模型（例如，低速气体和液体，例如水）。 直接手动求解这些方程并不容易，但很高兴我们可以获得简单问题的“精确”解。 使用这些方程来解决感兴趣的问题，比如流过机翼的空气或流过某个系统的水，需要我们对这些方程进行数值求解。

建立一个数值方案

为了使用计算机解决更复杂的问题，需要一种数值求解不可压缩方程的方法。 在数值上求解偏微分方程甚至微分方程并非易事。 然而，我们在本指南中的方程式有一个特殊的挑战（惊喜！）。 也就是说，我们需要求解动量方程，同时保持解散度不受连续性的要求。 通过像 Runge-Kutta 方法这样的简单时间积分变得困难，因为连续性方程中没有时间导数。

没有正确的，甚至是最好的求解方程的方法，但有许多可行的选择。 几十年来，已经发现了几种解决该问题的方法，例如根据涡度和流函数重新制定公式，引入人工可压缩性和算子分裂。 Chorin (1969)，然后是 Kim 和 Moin (1984, 1990)，制定了一种非常成功和流行的分步法，这将使我们能够在直接而不是隐式求解压力场的同时对方程进行积分。 分步法是一种通过分裂算子来逼近方程的一般方法，在这种情况下，沿着压力分裂。 该方法相对简单但稳健，促使其选择此处。

首先，我们需要在时间上对方程进行数值区分，以便我们可以从一个时间点跨步到下一个时间点。 继 Kim 和 Moin (1984) 之后，我们将对对流项使用二阶显式 Adams-Bashforth 方法，对粘性项使用二阶隐式 Crank-Nicholson 方法，对时间导数使用简单的有限差分，同时忽略压力梯度。 这些选择绝不是唯一可以做出的近似：它们的选择是通过控制方案的数值行为来构建方案的艺术的一部分。

中间速度现在可以积分，但是，它忽略了压力的贡献并且现在是发散的（不可压缩性要求它是无发散的）。 需要运算符的其余部分将我们带到下一个时间步。

在哪里是我们需要找到的一些标量，它会导致发散的自由速度。 我们可以找通过采取修正步骤的分歧，

其中第一项根据连续性的要求为零，产生标量场的泊松方程，该方程将在下一个时间步提供螺线管（无发散）速度。

正如 Kim 和 Moin (1984) 所示， 不完全是算子分裂造成的压力，但可以通过

在本教程的这一点上，我们做得很好，我们已经对控制方程进行了时间离散化，以便我们可以整合它们。 我们现在需要对算子进行空间区分。 有许多方法可以实现这一点，例如，有限元法、有限体积法和有限差分法。 在 Kim 和 Moin (1984) 的原创作品中，他们采用了有限差分法。 该方法因其相对简单和计算效率而具有优势，但由于需要结构化网格，因此会受到复杂几何形状的影响。

有限元法 (FEM) 因其通用性而成为一种方便的选择，并且有一些非常好的开源项目可以帮助其使用。 特别是它可以处理一维、二维和三维的真实几何图形，可以针对机器集群上的非常大的问题进行扩展，并且相对容易用于高阶元素。 通常，该方法是这三种方法中较慢的一种，但它会给我们最大的解决问题的里程，所以我们将在这里使用它。

即使在实施 FEM 时，也有很多选择。 在这里，我们将使用 Galerkin FEM。 在这样做时，我们通过将每个方程乘以一个测试函数来将方程转换为加权残差形式对于向量和对于标量场，并在域上积分. 然后我们使用斯托克定理或发散定理对任何高阶导数进行部分积分。 然后我们提出变分问题，产生我们想要的 CFD 方案。

我们现在有一个很好的数学方案，以“方便”的形式实施，希望能对实现目标有一定的了解（大量数学，以及来自杰出研究人员的方法，我们几乎复制和调整）。