srcmini本文概述
从历史上看, 计算科学在很大程度上局限于研究科学家和博士候选人的领域。但是, 这些年来, 也许对于大型软件社区而言并不为人所知, 但我们科学计算的笨蛋一直在悄悄地编译可处理绝大多数繁重工作的协作式开源库。结果是, 现在比以往任何时候都更容易实现数学模型, 并且尽管该领域可能尚未准备好进行大规模消费, 但成功实现这一目标的门槛已大大降低。从头开始开发新的计算代码库是一项艰巨的任务, 通常需要数年才能完成, 但是这些开放源代码科学计算项目使得可以通过可访问的示例来运行, 从而相对较快地利用这些计算能力。
由于科学计算的目的是提供对自然界中存在的真实系统的科学见解, 因此该学科代表了使计算机软件方法变为现实的最前沿。为了使软件能够在很高的精度和分辨率下模拟现实世界, 必须调用复杂的微分数学, 这需要大学, 国家实验室或公司研发部门所不具备的知识。最重要的是, 当试图使用零和一的离散语言描述现实世界的连续, 无穷小结构时, 就会遇到重大的数字挑战。为了使算法既具有计算上的可控性, 又能提供有意义的结果, 必须进行详尽的数值转换以进行详尽的工作。换句话说, 科学计算是繁重的工作。
科学计算的开源工具
我个人特别喜欢FEniCS项目, 将其用于我的论文工作, 并将证明我在选择本教程的代码示例中的偏见。 (还有其他一些非常高质量的项目, 例如DUNE也可以使用。)
FEniCS被自我描述为”一个协作项目, 旨在开发用于自动化科学计算的创新概念和工具, 并特别关注通过有限元方法自动解决微分方程的问题。”这是一个功能强大的库, 用于解决大量问题和科学计算应用程序。它的贡献者包括Simula研究实验室, 剑桥大学, 芝加哥大学, 贝勒大学和KTH皇家理工学院, 他们在过去十年中共同将其构建为宝贵资源(请参阅FEniCS规范)。
令人惊讶的是FEniCS库为我们提供了巨大的保护。为了了解该主题涵盖的深度和广度, 可以查看其开源教科书, 其中第21章甚至比较了解决不可压缩流的各种有限元方案。
该项目在幕后为我们集成了一大套开源科学计算库, 这些库可能很有趣, 也可以直接使用。其中包括FEniCS项目调用的项目(无特定顺序)：
- PETSc：一套数据结构和例程, 用于通过偏微分方程建模的科学应用程序的可扩展(并行)解决方案。
- Trilinos项目：一组强大的算法和技术, 用于解决线性和非线性方程式, 是根据Sandia国家实验室的工作开发的。
- uBLAS：”一个C ++模板类库, 为密集, 压缩和稀疏矩阵提供了BLAS 1、2、3级功能, 并为线性代数提供了许多数值算法。”
- GMP：一个用于任意精度算术的免费库, 可对有符号整数, 有理数和浮点数进行运算。
- UMFPACK：使用非对称MultiFrontal方法求解非对称稀疏线性系统Ax = b的一组例程。
- 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方程式, 事不宜迟：
质量守恒：
保持动量：
节约能源：
偏应力为：
尽管非常笼统, 它控制着物理世界中的大多数流体流, 但直接使用它们并没有多大用处。这些方程式的已知精确解相对较少, 并且对于任何能够解决存在性和平滑性问题的人, 这里都将获得100万美元的千年奖。重要的是, 我们通过进行一系列降低复杂性的假设(这是古典物理学中最难的一些方程式)来为模型的开发提供起点。
为了使事情”简单”, 我们将使用我们的领域特定知识对流体进行不可压缩的假设, 并假设温度恒定, 这样就不需要(成为热方程的)能量守恒方程(解耦)。现在, 我们有两个方程, 仍然是PDE, 但是在解决大量实际流体问题的同时, 方程式却明显简化了。
连续性方程
动量方程
此时, 我们现在有了一个很好的数学模型, 用于不可压缩的流体流动(例如低速气体和液体, 例如水)。直接手工求解这些方程式并不容易, 但是它可以为简单的问题获得”精确”的解决方案, 这是一个很好的选择。使用这些方程式来解决感兴趣的问题, 例如流过机翼的空气或流过某些系统的水, 要求我们用数值方法求解这些方程式。
建立数值方案
为了使用计算机解决更复杂的问题, 需要一种数值求解我们不可压缩方程的方法。在数值上求解偏微分方程甚至微分方程并非易事。但是, 本指南中的方程式有一个特殊的挑战(惊奇！)。就是说, 我们需要求解动量方程, 同时保持解的连续性, 从而使解无散度。由于连续性方程式内部没有时间导数, 因此很难通过Runge-Kutta方法进行简单的时间积分。
没有正确的, 甚至最好的方法来求解方程式, 但是有许多可行的选择。在过去的几十年中, 已经找到了解决该问题的几种方法, 例如在涡度和流函数方面重新制定公式, 引入人为可压缩性以及操作符拆分。 Chorin(1969), 然后是Kim和Moin(1984, 1990), 提出了一种非常成功且流行的分数阶跃方法, 该方法将使我们能够在直接求解压力场的同时(而不是隐式地)求解方程组。分数步长法是一种通用的方法, 可通过拆分等式(在这种情况下是沿压力拆分)来近似方程。该方法相对简单且健壮, 因此激发了它的选择。
首先, 我们需要在时间上对这些方程进行数值区分, 以便我们可以从一个时间点移至下一个时间点。在遵循Kim和Moin(1984)的基础上, 我们将对流项使用二阶显式Adams-Bashforth方法, 对粘性项使用二阶隐式Crank-Nicholson方法, 对时间导数进行简单的有限差分, 而忽略了压力梯度。这些选择绝不是唯一可以近似的方法：它们的选择是通过控制方案的数值行为来构建方案的技术的一部分。
中间速度现在可以积分, 但是它忽略了压力的作用, 现在发散了(不可压缩性要求它没有发散)。需要其余的运算符才能将我们带入下一个步骤。
其中
是我们需要找到的标量, 导致自由速度发散。我们可以找
通过采取纠正步骤的分歧,
其中第一项为0(如连续性所要求), 则产生标量场的泊松方程, 该泊松方程将在下一个时间步提供螺线管速度(发散自由)。
正如Kim和Moin(1984)所表明的,
不完全是由于操作员分裂而产生的压力, 但可以通过
.
至此, 在本教程中我们做得还不错, 我们暂时区分了控制方程, 以便我们可以对它们进行积分。现在, 我们需要在空间上区分运营商。存在许多可以实现此目的的方法, 例如, 有限元法, 有限体积法和有限差分法。在Kim和Moin(1984)的原始工作中, 他们采用了有限差分法。该方法因其相对简单和计算效率而具有优势, 但由于需要结构化的网格而使其几何形状复杂。
有限元方法(FEM)具有通用性, 是一种方便的选择, 并且有一些非常好的开源项目有助于其使用。特别是, 它可以处理一维, 二维和三维的实际几何形状, 可以针对机器集群上的非常大的问题进行缩放, 并且对于高阶元素相对易于使用。通常, 该方法是这三种方法中较慢的一种, 但是它将为我们提供最大的解决问题的机会, 因此我们将在这里使用它。
即使实施FEM, 也有许多选择。在这里, 我们将使用Galerkin FEM。为此, 我们将方程乘以测试函数, 从而将方程转换为加权残差形式
对于矢量和
标量字段, 并在域上进行积分
。然后, 我们使用斯托克定理或散度定理对任何高阶导数进行部分积分。然后, 我们提出变分问题, 得出我们所需的CFD方案。
现在, 我们以”方便”的形式实施了一个不错的数学方案, 希望能对到达那里需要些什么有所了解(很多数学, 以及出色的研究人员提供的很多方法, 我们都在进行复制和调整)。
评论前必须登录！
注册