Historically, 计算科学在很大程度上局限于研究科学家和博士候选人的领域. 然而，随着时间的推移——也许不为更大的软件社区所知—— us 科学计算书呆子 一直在悄悄地编译协作开源库，以处理绝大多数繁重的工作. 其结果是，现在实施起来比以往任何时候都容易 数学模型, 虽然这个领域可能还没有为大规模消费做好准备, 成功实施的门槛已经大大降低. 从头开始开发新的计算代码库是一项艰巨的任务, 通常以年为单位, 但是这些开放源码的科学计算项目使得使用可访问的示例来相对快速地利用这些计算能力成为可能.

因为科学计算的目的是提供对自然界中存在的真实系统的科学洞察, 该学科代表了使计算机软件接近现实的前沿. 为了使软件能够以非常高的精度和分辨率模拟现实世界, 必须调用复杂的微分数学, 需要在大学围墙之外很少能找到的知识, national labs, or corporate R&D departments. 最重要的是，当 试图描述现实世界中连续的，无限小的结构 使用0和1的离散语言. 为了使算法在计算上易于处理，需要进行详尽的数值变换, 并产生有意义的结果. 换句话说，科学计算任务繁重.

科学计算的开源工具

我个人特别喜欢 FEniCS project, 用它来写论文, 并将演示我在本教程的代码示例中选择它的偏见. (还有其他非常高质量的项目，如 DUNE 哪一个也可以用.)

FEniCS自我描述为“为自动化科学计算开发创新概念和工具的协作项目”, 特别侧重于用有限元方法自动求解微分方程.“这是一个强大的库，可以解决大量的问题和科学计算应用. 它的贡献者包括Simula研究实验室, 剑桥大学, 芝加哥大学, Baylor University, 和皇家理工学院, 在过去的十年中，他们共同将其建设成为一种宝贵的资源(参见FEniCS代码).

令人惊讶的是，FEniCS库为我们省去了多少工作. 要了解该项目涵盖的主题的惊人深度和广度，您可以 查看他们的开源教科书，其中第21章甚至比较了求解不可压缩流动的各种有限元格式.

在幕后，这个项目为我们集成了大量的开源科学计算库, 哪个可能是感兴趣的或直接使用的. 这些项目包括FEniCS项目提出的项目，排名不分先后:

• PETSc一套数据结构和例程，用于由偏微分方程建模的科学应用的可扩展(并行)解决方案.
• Trilinos Project一套用于求解线性和非线性方程的鲁棒算法和技术, 这是桑迪亚国家实验室的研究成果.
• uBLAS一个c++模板类库，提供BLAS级别1, 2, 3 .功能为密集, 填充矩阵和稀疏矩阵以及线性代数的许多数值算法.”
• GMP任意精度算术的自由库, 对有符号整数的运算, rational numbers, 还有浮点数.
• UMFPACK一组求解非对称稀疏线性系统的例程, Ax=b, 使用非对称多额方法.
• ParMETIS:一个基于mpi的并行库，实现了各种算法来划分非结构化图, meshes, 以及计算稀疏矩阵的减填充排序.
• NumPy:用Python进行科学计算的基本包.
• CGAL高效可靠的几何算法，以c++库的形式呈现.
• SCOTCH:用于顺序和并行图分区的软件包和库, 静态映射和集群, 顺序网格和超图分区, 以及顺序和并行稀疏矩阵块排序.
• MPI一种标准化的、便携的消息传递系统，由学术界和工业界的一组研究人员设计，可在各种并行计算机上运行.
• VTK: An open-source, 免费的三维计算机绘图软件系统, 图像处理和可视化.
• SLEPc:在并行计算机上求解大规模稀疏特征值问题的软件库.

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

同样有趣的是，除了科学计算之外，还有许多应用程序可以利用这些项目, 包括金融建模，图像处理，优化问题，甚至视频游戏. 这是可能的, for example, 创建一个视频游戏，使用这些开源算法和技术来解决二维流体流动问题, 比如玩家会与之互动的海洋/河流(游戏邦注:也许玩家会尝试着驾驶一艘船穿越不同的风和水流)。.

示例应用程序:利用开放源码进行科学计算

在这里，我将通过展示一个基本的计算流体动力学方案是如何在一个开源库中开发和实现的，来尝试给一个开发数值模型的感觉 FEniCS project. FEnICS提供Python和c++的api. 对于本例，我们将使用他们的Python API.

我们将讨论一些相当技术性的内容, 但我们的目标只是让大家体验一下开发这样的科学计算代码需要做些什么, 以及今天的开源工具为我们做了多少跑腿的工作. 在这个过程中，希望我们能帮助揭开科学计算这个复杂世界的神秘面纱. (Note that an Appendix 是否提供了所有的数学和科学基础的细节给那些对细节感兴趣的人.)

DISCLAIMER: 对于那些在科学计算软件和应用方面很少或没有背景的读者, 这个例子的某些部分可能会让你有这样的感觉:

如果是这样，不要绝望. 这里的主要收获是现有的开源项目可以极大地简化这些任务的程度.

记住这一点，让我们从FEnICS开始 不可压缩纳维-斯托克斯的演示. 本演示模拟了不可压缩流体流过l型弯道的压力和速度, 比如水管.

在演示页面上的描述链接给出了一个很好的, 运行代码的必要步骤的简明设置, 我鼓励大家快速浏览一下其中的内容. 综上所述，演示将求解通过弯道的速度和压力 不可压缩流动方程. 该演示运行了一个简短的流体随时间流动的模拟，并将结果动画化. 这是通过在管道中设置代表空间的网格，并使用 有限元法 对网格上各点的速度和压力进行数值求解. 然后我们通过时间进行迭代, 不断更新速度和压力场, 还是用我们手头的方程.

该演示按原样运行良好，但我们将对其稍加修改. The demo uses the Chorin split，但我们将使用稍微不同的方法，灵感来自 Kim and Moin 相反，我们希望它更稳定. 这只需要我们改变用来近似对流项和粘性项的方程, 但要做到这一点，我们需要存储前一个时间步长的速度场, 在更新方程中再加上两项, 哪一种会利用之前的信息得到更精确的数值近似.

我们来做这个改变. First, we add a new Function object to the setup. 这是一个表示抽象数学函数(如矢量或标量场)的对象. We’ll call it un1 它会存储之前的速度场

在函数空间上 V.

...
# Create functions	(三个不同的向量场和一个标量场)
un1 = Function(V)	#我们添加的前一个时间步长的速度场
u0 = Function(V)	#当前速度场
u1 = Function(V)	#下一个速度场(正在解的是什么)
p1 = Function(Q)	#下一个压力场(正在解决的问题)
...

Next, 我们需要改变在模拟的每个步骤中更新“暂定速度”的方式. 这个场表示忽略压力时下一个时间步长的近似速度(此时压力尚不可知). 这是我们用最近的Kim和Moin分步法取代Chorin分割法的地方. 换句话说，我们将更改字段的表达式 F1:

Replace:

#暂定速度场(对下一个速度场的第一个预测)
#为Chorin风格分割
# F1 =速度场变化量+
#对流项+
#扩散项
#身体力项

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

With:

#暂定速度场(对下一个速度场的第一个预测)
# Kim和Moin风格的分裂
# F1 =速度场变化量+
#对流项+
#扩散项
#身体力项

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，在每个迭代步骤结束时

...
移动到下一个时间步骤
un1.assign(u0)		#复制当前的速度场到之前的速度场
u0.assign(u1)		#复制下一个速度场到当前速度场
...

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

这个演示程序解决了不可压缩的Navier-Stokes方程
用Kim和Moin的分步法在l形域上进行了研究."""

# Begin demo

从dolfin导入*

#只从根进程并行打印日志消息
parameters["std_out_all_processes"] = False;

#从文件中加载网格
mesh = Mesh("lshape.xml.gz")

#定义函数空间(P2-P1)
V = VectorFunctionSpace(mesh， "Lagrange"， 2)
Q = FunctionSpace(mesh， "Lagrange"， 1)

定义trial和test函数
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

#设置参数值
dt = 0.01
T = 3
nu = 0.01

定义随时间变化的压力边界条件
p_in =表达式("sin(3.0*t)", t=0.0)

#定义边界条件
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 =[流入，流出]

# Create functions
un1 = Function(V)
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)

#定义系数
k = Constant(dt)
f = Constant((0,0))

#暂定速度场(对下一个速度场的第一个预测)
# Kim和Moin风格的分裂
# F1 =速度场变化量+
#对流项+
#扩散项
#身体力项

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)

如果可用，使用amg预调节器
Prec = "amg" if has_krylov_solver_preconditioning ("amg") else "default"

#创建存储解决方案的文件
ufile = File("results/velocity。.pvd")
pfile = File("结果/压力.pvd")

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

    #更新压力边界条件
    p_in.t = t

    计算暂定速度步长
    开始(“计算暂定速度”)
    b1 = assemble(L1)
    [bc.在bcu中申请(A1, b1) bc]
    solve(A1, u1.Vector ()， b1， "gmres"， "default")
    end()

    #压力校正
    开始(“计算压力校正”)
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    solve(A2, p1.Vector ()， b2， "cg"， prec)
    end()

    #速度校正
    begin(“计算速度校正”)
    b3 = assemble(L3)
    [bc.申请(A3, b3)的bc在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

    移动到下一个时间步骤
    un1.assign(u0)
    u0.assign(u1)
    t += dt
    print "t =", t

# Hold plot
interactive()

运行程序会显示弯头周围的流量. 自己运行科学计算代码，看看它的动画效果! 最后一帧的屏幕如下所示.

模拟结束时弯道内的相对压力, 按大小缩放和着色(无量纲值):

模拟结束时弯道中的相对速度，以矢量符号的形式按大小缩放和着色(无量纲值).

我们所做的是一个现有的演示, 哪一种方案恰好很容易实现与我们的方案非常相似, 然后利用前一个时间步长的信息对它进行了修改，使其具有更好的近似值.

在这一点上，您可能认为这是一个微不足道的编辑. 是的，这就是问题的关键所在. 这个开源的科学计算项目允许我们通过更改四行代码来快速实现修改后的数值模型. 在大型研究代码中，这样的改变可能需要几个月的时间.

该项目还有许多其他的演示，可以作为一个起点. 甚至有很多 开源应用程序 构建在项目上实现各种模型.

Conclusion

科学计算及其应用确实很复杂. 这是无法回避的. 但在许多领域，这是越来越正确的, 不断增长的可用开源工具和项目可以显著简化原本极其复杂和乏味的编程任务. 也许科学计算变得足够容易使用的时候已经不远了，它可以很容易地在研究社区之外被使用.

---

附录:科学和数学基础

有兴趣的同学, 以下是我们上面的计算流体动力学指南的技术基础. 以下内容将作为一个非常有用和简洁的主题总结，这些主题通常涵盖在十几个研究生水平的课程中. 研究生和数学类型有兴趣深入了解这个主题可能会发现这个材料相当吸引人.

Fluid Mechanics

一般来说，“建模”是用一系列近似解实际系统的过程. The model, 通常会涉及不适合计算机实现的连续方程, 因此必须用数值方法进一步逼近.

对于流体力学, 让我们从基本方程开始, Navier-Stokes方程, 并用它来开发CFD方案.

The n - s方程 are a series of 偏微分方程(PDEs) 它们很好地描述了流体的流动，因此是我们的起点. 它们可以从质量、动量和能量守恒定律推导出来 雷诺输运定理 and applying Gauss’s theorem 并引用斯托克假说. 这些方程需要连续介质假设, 假设我们有足够的流体粒子来给出统计性质，比如温度, 密度和速度含义. 此外，表面应力张量与应变速率张量之间存在线性关系, 应力张量的对称性, 各向同性流体的假设是必要的. 了解我们在开发过程中做出和继承的假设是很重要的，这样我们就可以评估结果代码中的适用性. Navier-Stokes方程 Einstein notation，废话少说:

质量守恒:

动量守恒:

能量守恒:

其中，偏应力为:

While very general, 控制着物质世界中大多数流体的流动, 它们没有多少直接用途. 已知的精确解相对较少，a 千禧年奖金100万美元 有没有人能解决 存在性和平滑性问题. 重要的是我们有了一个发展模型的起点通过做出一系列的假设来减少复杂性(它们是经典物理学中最难的一些方程).

为了使事情“简单”，我们将使用我们的领域特定知识来对流体进行不可压缩假设, 假设温度恒定，能量守恒方程, 哪个变成了热方程, 不需要(解耦). 现在我们有两个方程, still PDEs, 但在解决大量实际流体问题的同时，要简单得多.

连续性方程

Momentum equations

至此，我们有了不可压缩流体流动(低速气体和水等液体)的数学模型, for example). 直接用手解这些方程并不容易, 但它的好处在于我们可以得到简单问题的“精确”解. 用这些方程来解决感兴趣的问题, 比方说气流流过机翼, 或者水流经某个系统, 要求我们用数值方法解这些方程.

建立数值格式

为了用计算机解决更复杂的问题, 需要一种数值解不可压缩方程的方法. 用数值方法求解偏微分方程，甚至微分方程，并不简单. 然而，我们在本指南中的方程式有一个特殊的挑战(惊喜)!). That is, 我们需要在解散度自由的情况下解动量方程, 作为连续性的要求. 一个简单的时间积分 Runge-Kutta method 是不是因为连续性方程里面没有时间导数而变得很困难.

没有正确的, or even best, 解方程的方法, 但是有很多可行的选择. 几十年来，人们找到了几种方法来解决这个问题, 例如根据涡度和流函数重新表述, 引入人工压缩性, 算子分裂. Chorin (1969), 然后是Kim和Moin (1984), 1990), 提出了一种非常成功和流行的分步法，使我们可以在直接求解压力场的同时对方程进行积分, 而不是含蓄地. 分步法是一种通过分割算子逼近方程的一般方法, 在这种情况下，沿着压力分裂. 这种方法相对简单，但却很健壮，这促使我们在这里选择它.

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

中间速度现在可以积分了, however, 它忽略了压力的作用，现在是发散的(不可压缩性要求它无发散). 运算符的余数需要带我们到下一个时间步长.

where

我们需要找一个标量来得到发散自由速度吗. We can find 通过修正步的散度，

第一项是0，这是连续性的要求，得到a Poisson equation 对于一个标量场，它将在下一个时间步提供一个螺线形(发散自由)速度.

正如《欧博体育app下载》(1984)所展示的那样，

压力不完全是由于操作者的分裂，而是可以由它求得的吗

.

在教程的这一点上，我们做得相当不错, 我们对控制方程进行了时间判别，以便对它们进行积分. 我们现在需要在空间上区分算子. 我们有很多方法可以做到这一点 有限元法, the 有限体积法, and the 有限差分法, for example. 在Kim和Moin(1984)的原著中，他们继续使用有限差分法. 该方法具有相对简单和计算效率高的优点, 但是对于复杂的几何形状来说，它需要 structured mesh.

有限元法(FEM)是一种方便的选择，因为它具有通用性，并且有一些非常好的开源项目帮助它的使用. 特别地，它可以一次性处理真实的几何图形, 二维和三维, 适用于机器集群上非常大的问题, 对于高阶元素来说相对容易使用. 通常这种方法是三种方法中较慢的一种, 但是它会给我们最大的跨越问题的里程, 我们在这里用它.

即使在实施FEM时，也有许多选择. 这里我们将使用 Galerkin FEM. 在这样做时，我们通过将每个方程乘以一个测试函数来将方程转换为加权残差形式

对于向量和 对于标量场，在定义域上积分 . 对任意高阶导数进行偏积分 Stoke’s theorem or the Divergence theorem. 然后我们提出变分问题，得到我们想要的CFD方案.

现在我们有了一个“方便”的数学方案, 希望你能对实现这一目标所需要的东西有所了解(大量的数学), 以及来自杰出研究人员的方法，我们基本上是复制和调整).