UFL：有限元形式语言》算法(二)&实现的问题【翻译】

UFL：有限元形式语言》算法(二)&实现的问题【翻译】

2021-02-04
| 科学计算 | | FEniCS , 有限元法 , 偏微分方程 | Comment 评论

【章节目录】

17·8 算法（续前）

17·8·4 重要的变换

有很多方法可以操纵表达式的表示。 在这里，我们描述了三个特别重要的变换。 注意，这里每一个算法都删除了一些抽象，因此可能删除了一些分析或优化的机会。 为了展示其效果，下面将每种变换应用到下面的表达式

\[ a = \mathrm{grad}(f u) \cdot \mathrm{grad}\ v \tag{17.77} \]

在本节的最后，给出了一些示例代码来演示更多的表示细节。

UFL中的某些算符称为“复合”算符，这意味着它们可以由其他更基本的算符表示。 尝试定义一个表达式a = dot(grad(f*u), grad(v))，然后打印repr(a)。 如您所见，a的表示是Dot(Grad(Product(f, u)), Grad(v))，其中一些详细信息代替了f，u和v。 通过直接使用高级别类型Grad（而不是更多低级别类型）来表示梯度，输入表达式在表示中更容易识别，并且将表达式呈现为比如LATEX格式可以显示终端用户所写的原始复合算符。 但是，由于许多算法必须为每种算符类型实现操作，因此函数expand_compounds可以将所有“ 复合”类型的表达式节点替换为使用基本类型的等效表达式。 当将此操作应用于来自用户的输入形式时，UFL和形式编译器中的算法仍可以纯粹用更基本的算符来编写。 展开方程式（17.77）的复合表达式可得到表达式

\[ a_c = \sum_i\frac{\partial v}{\partial x_i}\frac{\partial (uf)}{\partial x_i} \tag{17.78} \]

另一个重要的变换是expand_derivatives，递归地且对所有类型的导数，应用对表达式的自动微分。 最终结果是大多数导数都能被计算，并且表达式树中剩余的仅有导数算符类型是应用到终端表达式的。 该算法的前提是已应用过expand_compounds。 根据公式（17.78）对ac展开求导可得出

\[ a_d = \sum_i\frac{\partial v}{\partial x_i}\left(u \frac{\partial f}{\partial x_i}+ f \frac{\partial u}{\partial x_i}\right) \tag{17.79} \]

索引记号和IndexSum表达式节点类型在某种程度上使表达式树的解释复杂化，尤其是在具有嵌套索引求和的表达式中。 由于具有自由索引的表达式将具有多个值，因此每个表达式对象不仅代表一个值，而且代表一组值。 变换expand_indices便派上用场了。 该算法的前提是已应用过expand_compounds和expand_derivatives。 该算法的后置条件是表达式中没有剩余的自由索引。 最终将公式（17.79）展开索引，得到

\[ a_i = \frac{\partial v}{\partial x_0}\left(u\frac{\partial f}{\partial x_0}+ f \frac{\partial u}{\partial x_0}\right) +\frac{\partial v}{\partial x_1}\left(u\frac{\partial f}{\partial x_1}+ f \frac{\partial u}{\partial x_1}\right) \tag{17.80} \]

我们从等式（17.77）中的更高层次的概念梯度和点积开始，最后仅以形式参数的标量加法，乘法和偏导数结束。 形式编译器通常以ad或ai开头，插入参数导数值，并在最终生成代码之前应用一些其他变换。

一些示例代码可以帮助您理解这些算法在表达式的表示级别上的作用。 由于此代码的打印输出有点长，因此下面仅输出要重复的关键方面。 将此代码复制到python文件或在python解释器中运行以查看完整的输出。

# Python code

from ufl import *
V = FiniteElement("Lagrange", triangle, 1)
u = TestFunction(V)
v = TrialFunction(V)
f = Coefficient(V)

# Note no *dx! This is an expression, not a form.
a = dot(grad(f*u), grad(v))

from ufl.algorithms import *
ac = expand_compounds(a)
ad = expand_derivatives(ac)
ai = expand_indices(ad)
print("\na: ", str(a), "\n", tree_format(a))
print("\nac:", str(ac), "\n", tree_format(ac))
print("\nad:", str(ad), "\n", tree_format(ad))
print("\nai:", str(ai), "\n", tree_format(ai))

显示为 \(a\) 的打印输出是（有限元对象的详细信息，已切除部分信息，以缩短行数）：

# Output
a: (grad(v_{-2} * w_0)) . (grad(v_{-1}))
Dot
(
    Grad
        Product
        (
            Argument(FiniteElement(...), -2)
            Coefficient(FiniteElement(...), 0)
        )
    Grad
        Argument(FiniteElement(...), -1)
)

标为-1和-2的参数分别引用v和u。

在 \(a_c\) 中，Dot乘法已展开为具有两个Indexed操作数的Product的IndexSum：

# Output
IndexSum
(
    Product
    (
        Indexed
        (
            ...
            MultiIndex((Index(10),), {Index(10): 2})
        )
        Indexed
        (
            ...
            MultiIndex((Index(10),), {Index(10): 2})
        )
    )
    MultiIndex((Index(10),), {Index(10): 2})
)

看起来有点复杂的表达式MultiIndex((Index(10),), {Index(10): 2})可以简单地理解为“绑定到第2维轴的名为 \(i_{10}\) 的索引”。

放大到上面的…行之一，将grad(f u)的表示变换为更多基本表达式后，仍必须保持矢量形状，这就是为什么SpatialDerivative对象包装在ComponentTensor对象中的原因： 【译者注：较新版本显示的信息和这里不一样，具体以实际打印信息为准，后面不再作此说明】

# Output
ComponentTensor
(
    SpatialDerivative
    (
        Product
        (
            u
            f
        )
        MultiIndex((Index(8),), {Index(8): 2})
    )
    MultiIndex((Index(8),), {Index(8): 2})
)

在展开表达式的算法中，出现的常见模式：

# Output
Indexed
(
    ComponentTensor
    (
        ...
        MultiIndex((Index(8),), {Index(8): 2})
    )
    MultiIndex((Index(10),), {Index(10): 2})
)

此模式用作索引对象的重新标记，将. . .所再内部的 \(i_8\) 重命名为外部的 \(i_{10}\) 。 当查看 \(a_d\) 的打印结果时，链式规则 \(((f u)′ = u f ′ + f u′)\) 的结果可以看作是两个Product对象的Sum。

# Output
Sum
(
    Product
    (
        u
        SpatialDerivative
        (
            f
            MultiIndex((Index(8),), {Index(8): 2})
        )
    )
    Product
    (
        f
        SpatialDerivative
        (
            u
            MultiIndex((Index(8),), {Index(8): 2})
        )
    )
)

最终，在 \(a_i\) （此处未显示）中进行索引展开之后，没有剩余自由的Index对象，反而在 \(a_i\) 的打印中看到很多FixedIndex对象。 如果您想很好地了解此处显示的三个变换，强烈建议您仔细阅读上面示例代码的全部输出。

17·8·5 表达式求值

即使UFL表达式打算由形式编译器进行编译，将其直接求值得浮点值也很有用的。 特别是，这使得UFL的测试和调试更加容易，并且在单元测试中得到了广泛的使用。 要对UFL表达式求值，必须指定形式参数和几何量的值。 通过将带有坐标的元组传递给调用算符，可对只依赖于空间坐标的表达式求值。 可以直接复制到交互式Python会话的如下代码展示了这种语法：

# Python code
from ufl import *
cell = triangle
# 在较新的版本中，cell没有x这个成员，改用随后的代码
# x = cell.x
x = SpatialCoordinate(cell)
e = x[0] + x[1]
print(e((0.5, 0.7))) # prints 1.2

可以使用从终端表达式实例到值的映射字典来指定其他终端表达式。 此代码通过映射扩展了上面的代码：

# Python code
c = Constant(cell)
e = c*(x[0] + x[1])
print(e((0.5, 0.7), { c: 10 })) # prints 12.0

如果函数和基函数依赖于空间坐标，那么这个映射可以指定一个Python可调用函数，而不是字面常量。 这个可调用函数必须是将空间坐标作为输入并返回一个浮点值。 如果要映射的函数是向量函数，则可调用函数必须返回值的元组。 这些扩展可以在以下代码中看到：

# Python code
element = VectorElement("Lagrange", triangle, 1)
c = Constant(triangle)
f = Coefficient(element)
e = c*(f[0] + f[1])
def fh(x):
    return (x[0], x[1])
print(e((0.5, 0.7), { c: 10, f: fh })) # prints 12.0

为了使用表达式求值来验证导数计算是否正确，还可以指定形式参数的空间导数。 然后，可调用函数必须携带第二个参数，该参数是指定要微分的空间方向的整数元组。 最终的示例代码，对 \(g = x_0 x_1\) ，计算出 \(g^2 + g^2_{，0} + g^2_{，1}\) 。

# Python code
element = FiniteElement("Lagrange", triangle, 1)
g = Coefficient(element)
e = g**2 + g.dx(0)**2 + g.dx(1)**2
def gh(x, der=()):
    if der == (): return x[0]*x[1]
    if der == (0,): return x[1]
    if der == (1,): return x[0]
print(e((2, 3), { g: gh })) # prints 49

17·8·6 表达式查看

可以用各种方式格式化表达式以进行检查，这在调试时特别有用。 内置的Python字符串转换算符str(e)提供了紧凑的人类可读字符串。 如果在交互式Python会话中键入print(e)，则显示str(e)。 另一个Python内置字符串算符是repr(e)。 UFL所正确实现的repr，满足对任何表达式e，e == eval(repr(e))都成立。 字符串repr(e)反映表达式中使用的所有确切表示类型，因此对于调试很有用。 另一个格式化函数是tree_format(e)，它产生一个缩进的多行字符串，该字符串清楚地显示了表达式的树结构，而repr则相反，后者可能返回很长且难以阅读的字符串。 可以在手册中找到有关将表达式格式化为LATEX和点图可视化格式的信息。

17·9 实现的问题

17·9·1 Python作为领域专用语言的基础

将UFL作为一个嵌入在Python的语言， 最初的这种选择会影响本节中详细介绍的许多实现细节。 因此，这里有一些关于为什么Python适用于此，为什么不适用的说法。

Python提供了一种简单的语法，通常被认为类似于伪代码。这是一个领域专用语言的很好起点。 很好地支持面向对象和操作符重载，这是UFL设计的基础。 Python的函数式编程功能（例如生成器表达式）在算法和形式编译器的实现中很有用。 内置的数据结构列表，字典和集合在可伸缩算法的快速实现中起着核心作用。

Python中的运算符重载有一个问题，那就是比较运算符。 问题源于以下事实：__eq__或__cmp__被用于内置数据结构字典，并将其设置为比较键，这意味着a == b必须返回布尔值，即使Expr也是如此。 结果是__eq__不能被重载，用来返回某些Expr类型的表示（比如Equals(a, b)）供形式编译器稍后处理。 另一个问题是and和or也不能被重载，因此不能在条件表达式中使用。 在Python中，这些设计选择是有充分的理由。 而这又是UFL中比较运算符设计有些不直观的原因。

17·9·2 确保形式签名的唯一性

形式编译器需要计算每个形式的唯一签名，以便在缓存系统中使用，避免重新编译。 定义签名的便捷方法是使用repr(form)，因为在Python中的此定义是eval(repr(form)) == form。 因此，对所有Expr子类都实现了__repr__。

有些形式在数学上是等效的，即使它们的表示不完全相同。 UFL的表达式没有使用真正的规范形式，而是采取了一些措施来确保识别这些平凡的等价形式。

Expr类层次结构中的某些类型（Counted的子类）具有全局计数器，用于标识创建它们的顺序。 形式参数（Argument和Coefficient）都使用了此计数器来标识它们在形式参数列表中的相对顺序。 其他要计数类型有Index和Label，它们仅将计数器用作唯一标识符。 实现了对所有Counted类型进行重新编号的算法，以使所有计数均从0开始。

另外，某些算符类型（例如Sum和Product）要维护操作数的排序列表，使得a + b和b + a都表示为Sum(a, b)。 此操作数排序有意与索引编号无关，因为那样会不稳定。 这种不稳定的原因是索引重编号算法的结果取决于操作数的顺序。 操作数排序和重编号相结合，可确保相等的形式的签名保持不变。 请注意，形式的表示以及签名可能会随UFL版本的变化而变化。 下面的行打印了应用了expand_derivatives并重新编号的形式签名。

# Python code
print(repr(preprocess(myform).preprocessed_form))

17·9·3 效率考量

在Python中编写UFL，我们显然没有将最高性能作为第一要务。 如果形式编译过程可以融合到应用程序构建过程中，那么性能就足够了。 但是，我们确实关心可扩展性能以有效地处理复杂的方程，因此，我们关心的是所使用算法的渐近复杂性。

要用Python编写清晰有效的算法，正确使用内置数据结构非常重要。 这些数据结构特别包括列表，字典和集合。 CPython（van Rossum等人），Python的参考实现，将数据结构列表实现为数组，这意味着append和pop以及随机读写访问都是O(1)操作。 但是，随机插入为O(n)。 字典和集合都实现为哈希映射，后者仅不具有与键关联的值而已。 在哈希映射中，只要键类型有效地实现了__hash__和__eq__，随机读写，插入和删除都是O(1)操作。 字典数据结构被Python语言广泛使用，因此要特别注意使其高效（Kuchling，2007年）。 因此，要享受高效地利用这些容器，所有Expr子类都必须有效地实现这两个特殊函数。 如此考虑对于使UFL的有效实现非常重要。

17·10 结论和未来方向

更多的其它功能可以引入到UFL。 添加哪些功能将取决于FEniCS用户和开发人员的需求。 某些功能可以单独在UFL中实现，但是大多数功能将需要更新FEniCS项目的其他部分。 因此，UFL的未来方向与整个FEniCS项目的发展紧密相关。

在UFL中，很容易对有限元的声明进行改进。 附加的复杂性主要是在形式编译器中。 当前的建议包括时空单元和时间导数。 具有不均匀胞元类型的其他几何映射和有限元空间也是有可能扩展。

可以添加其他算符以使语言更具表现力。 一些算符很容易添加，因为它们的实现只影响一小部分代码。 很容易添加更多的使用基本算符表示的复合算符。 只要知道其导数，其他特殊函数也很容易添加。 其他功能可能需要更全面的设计考量，例如对复数的支持，这会影响大部分代码。

用户友好的记号和对快速开发的支持是UFL设计的核心价值。 使用接近数学抽象符号可以更轻松地表达特定的想法，从而可以减少用户代码中错误的可能性。 但是，元编程和代码生成的概念增加了另一层抽象，这可能使终端用户对框架的理解更加困难。 因此，每处都进行良好的错误检查非常重要，尽可能地检测用户输入错误。 错误消息，文档和单元测试套件的改进将始终对您有所帮助，以避免新用户频繁重复出现错误和误解。

为了支持形式编译器项目，可以在UFL中包含用于更有效地生成更好代码的算法和实用程序。 这样的算法可能应该局限于诸如表达图的一般变换之类的算法，独立于形式编译器特定方法。 在这一领域，更多的关于自动微分算法替代方案的工作（Forth等，2004； Tadjouddine，2008）可能有用。

总而言之，UFL是FEniCS框架的核心组件，在此框架中，它提供了丰富的形式语言，自动微分功能以及高效形式编译器的构建块。 这些对于快速求解偏微分方程的应用程序很有用。 UFL提供了离散自动化的用户接口，这是FEniCS的核心功能，并将线性自动化添加到了框架中。 凭借这些功能，UFL使FEniCS距其总体目标“自动化数学建模”又近了一步。

17·11 致谢

这项工作得到了挪威研究委员会（grant 162730）和Simula研究实验室的支持。 我要感谢通过建议和测试帮助改善UFL的所有人，特别是Anders Logg，KristianØlgaard，Garth Wells，Harish Narayanan和Marie Rognes。 除了两位匿名评论， Kent-André Mardal和Marie Rognes都进行了严格的评论，极大地改善了本章。

【章节目录】