title: "UFL：有限元形式语言》形式算符&表达式的表示【翻译】"
date: 2021-02-02T12:56:05+08:00
draft: true
categories: ["科学计算"]
tags: ["FEniCS","有限元法","偏微分方程"]

【[[docs/fem/_index#第十七章 UFL：有限元形式语言|章节目录]]】

17·5 形式算符

一旦定义了一些形式，就有几种方法可以从中计算相关的形式。 上节的算符可用于定义表达式，本节中讨论的算符被用于形式，从而生成新的形式。 形式算符既可以使形式的定义更紧凑，又可以减少错误的可能，因为原始形式中的更改将自动传播到根据它所计算出的形式中。 这些形式算符可以任意组合； 给定一个半线性形式，只需要几条行即可计算出雅可比伴随的作用。 由于这些计算是在形式编译器处理之前完成的，因此在运行时没有任何开销。

17·5·1 形式的微分

形式算符derivative声明了一个形式关于离散函数系数（Coefficient）的导数。 例如，可以使用此功能结合Newton-Raphson方法，自动将您的非线性残差方程（线性形式）线性化。 它也可以被多次应用，这对于从凸函数导出线性系统很有用，以便找到使泛函最小化的函数。 对于非平凡方程，手工计算这些表达式可能很繁琐。 此功能可能对其他领域也是有用的，包括最优控制和反演理论，以及灵敏度分析。

就简单的形式，声明形式$L$关于系数函数$w$的导数：

# UFL code
a = derivative(L, w, u)

形式$a$取决于附加的基函数参数$u$，该参数必须与函数$w$位于相同的有限元空间中。 如果省略最后一个参数，则会创建一个新的基函数参数。

让我们逐步介绍一个示例，如何对泛函求导两次，进而导出一个线性系统。 在下文中，$V_h$是具有基底的有限元空间，$w$是$V_h$中的函数，而$f = f (w)$是我们要最小化的泛函。 从$f (w)$导出的是线性形式$F(w; v)$和双线性形式$J(w; u, v)$。

$$ V_h = \mathrm{span} {\phi_k} \tag{17.52} $$

$$ w(x) =\sum^{|V_h|}_{k=1}w_k \phi_k(x) \tag{17.53} $$

$$ f : V_h \to \mathbb{R} \tag{17.54} $$

$$ F(w; \phi_i) =\frac{\partial f (w)}{\partial w_i} \qquad i = 1, \dots , |V_h| \tag{17.55} $$

$$ J(w; \phi_j, \phi) =\frac{\partial F(w; \phi)}{\partial w_j} \qquad j = 1, \dots , |V_h|, \phi \in V_h \tag{17.56} $$

对于具体的泛函$f (w) = \int_\Omega \frac{1}{2} w^2 dx$，我们可以将其实现为

# UFL code
v = TestFunction(element)
u = TrialFunction(element)
w = Coefficient(element)
f = 0.5*w**2*dx
F = derivative(f, w, v)
J = derivative(F, w, u)

此代码声明两种形式F​和J。 线性形式F表示标准载荷矢量w*v*dx，双线性形式J表示质量矩阵u*v*dx。

也可以是关于有限元混合空间函数的系数的求导。 考虑从泛函导出的调和映射方程

$$ f (x, \lambda) = \int_\Omega \left(\mathrm{grad} \\ x : \mathrm{grad} \\ x + \lambda x \cdot x\right) dx \tag{17.57} $$

其中，$x$是向量有限元空间$V_h^d$中的函数，而$\lambda$是标量有限元空间$V_h$中的函数。 从等式17.57中的泛函导出的线性和双线性形式，有属于混合空间$V_h^d\times V_h$中的基函数参数。 可通过自动线性化来得到这些形式

# UFL code
Vx = VectorElement("Lagrange", triangle, 1)
Vy = FiniteElement("Lagrange", triangle, 1)
u = Coefficient(Vx*Vy)
x, y = split(u)
f = inner(grad(x), grad(x))*dx + y*dot(x,x)*dx
F = derivative(f, u)
J = derivative(F, u)

注意，该泛函通过两个子函数$x$和$y$来表示，而derivative的参数必须是单个混合函数$u$。 在此示例中，derivative的基函数参数被省略，因此自动由正确的函数空间中提供。

请注意，在计算形式的导数时，我们假设

$$ \frac{\partial }{\partial w_k} \int_\Omega I dx = \int_\Omega \frac{\partial}{\partial w_k} I dx \tag{17.58} $$

或更特别地，域$\Omega$与$w$无关。 同样，除$w$以外的任何系数都被认为与$w$无关。 此外，请注意，在此框架中对单元的选择没有任何限制，特别是支持任意混合单元。

17·5·2 伴随

另一个形式算子是双线性形式$a$的伴随$a^∗$，定义为$a^∗(v, u) = a(u, v)$，等效于组装成稀疏矩阵后再进行转置。 在UFL中，这可以简单地通过交换测试函数和试探函数的顺序来实现，并且可以使用形式算符adjoint。 （请注意，这不是伴随算符的最一般的定义）。 在各向异性扩散项上使用它的一个例子看起来像

# UFL code
V = VectorElement("Lagrange", cell, 1)
T = TensorElement("Lagrange", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
M = Coefficient(T)
a = M[i,j]*u[k].dx(j)*v[k].dx(i)*dx
astar = adjoint(a)

对应于（$u\in U$和$v\in V$）

$$ a(M; u, v) =\int_\Omega M_{ij} u_{k,j} v_{k,i} dx \tag{17.59} $$

$$ a^∗(M; v, u) =\int_\Omega M_{ij} u_{k,j} v_{k,i} dx = a(M; u, v) \tag{17.60} $$

如果我们需要使用求导来计算非对称双线性形式的伴随，那么这种自动转换特别有用，因为$a$的显式表达式并不在其间。 与derivative结合使用时，以下几种形式算符最有用。

17·5·3 替换函数

通过使用其他值替换终端对象，可以使用新定义的形式参数来对形式求值。 假定您已经定义了依赖于某些函数$f$和$g$的形式$L$。 然后，您可以通过将这些函数替换为其他函数或固定值，来得到特定形式，比如

$$ L(f , g; v) = \int_\Omega \frac{f^2}{2g}v\\ dx \tag{17.61} $$

$$ L_2(f , g; v) = L(g, 3; v) = \int_\Omega \frac{g^2}{6}v\\ dx \tag{17.62} $$

此功能通过replace实现了，如下所示：

# UFL code
V = FiniteElement("Lagrange", cell, 1)
v = TestFunction(V)
f = Coefficient(V)
g = Coefficient(V)
L = f**2 / (2*g)*v*dx
L2 = replace(L, { f: g, g: 3})
L3 = g**2 / 6*v*dx

L2和L3代表完全相同的形式。 由于它们仅取决于g，为这些形式生成的代码可以更高效。

17·5·4 作用

在某些应用中，不明确需要矩阵，而只需要矩阵对矢量的作用。 直接组装结果向量，先比组装稀疏矩阵然后执行矩阵-向量乘法要有效得多。 假设$a$是双线性形式，$w$是定义在与$a$中试探函数相同的有限元上的系数。 令$A$表示可以从$a$组装的稀疏矩阵。 然后，您可以通过定义代表双线性形式$a$作用于一个函数$w$的线性形式$L$， 进而直接将其组装成$A$对向量的作用。 简单地表示成L = action(a, w)，甚至缩写为L = a*w。

17·5·5 系统分割

如果您更喜欢将PDE的所有项都写在一侧，比如

$$ a(u, v) − L(v) = 0 \tag{17.63} $$

您可以同时声明线性项和双线性项，然后将方程拆分为$a$和$L$。 一个简单的例子是

# UFL code
V = FiniteElement("Lagrange", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Coefficient(V)
pde = u*v*dx - f*v*dx
a, L = system(pde)

这里的system用于将PDE分为双线性和线性部分。 或者，可以使用lhs和rhs分别获得这两个部分。 记住，线性部分的结果符号，对应于在公式（17.63）中将L移动到右侧。

17·5·6 计算函数的敏感性

如果您找到了方程（17.63）的解$u$，并且$u$依赖于某个标量常值$c$，则可以计算$u$关于$c$的变化灵敏度。 如果$u$被表示成代数线性系统$Ax = b$解的系数向量$x$，那么$\frac{\partial u}{\partial c}$的系数就是$\frac{\partial x}{\partial c}$。 将$\frac{\partial}{\partial c}$应用于$Ax = b$并使用链式规则，我们可以写成

$$ A \frac{\partial x}{\partial c}=\frac{\partial b}{\partial c}−\frac{\partial A}{\partial c} x \tag{17.64} $$

因此，可以通过求解相同的代数线性系统算出的$x$来发现$\frac{\partial x}{\partial c}$， 仅在右侧不同。 可以写出与等式（17.64）右边相对应的线性形式

# UFL code
u = Coefficient(element)
sL = diff(L, c) - action(diff(a, c), u) # 执行报错，我尚未做到原因

或者您可以使用等效的形式转换

# UFL code
sL = sensitivity_rhs(a, u, L, c) # 执行报错，我尚未做到原因

请注意，解$u$必须用Coefficient表示，而$a(u, v)中$的$u$用Argument表示。

17·6 表达式的表示

从高级的角度来看，UFL就是形式的定义。 每种形式都包含一个或多个标量被积表达式，但是形式表示在很大程度上与被积表达式的表示无关。 实际上，UFL实现的大多数复杂性都与表达式的表示，表现和操纵有关。 本章的其余部分将重点介绍表达式的表示以及在其上操作算法。 这些主题对UFL的普通用户几乎没有兴趣，而更多地针对开发人员和好奇的面向技术的用户。

为了在没有实现细节负担的情况下推理表达式算法，我们需要对表达式结构进行抽象表示。 UFL表达式是程序的表示，并且该表示法应使我们能够看到这种联系。 下面，我们将根据这种抽象符号来讨论表达式的属性，并与特定的实现细节相关。

17·6·1 表达式的结构

不依赖于其他表达式的最基本的表达式称为终端表达式。 将某些算符应用于一个或多个现有表达式会产生其他表达式。 考虑任意（非终端）表达式$z$。 该表达式依赖于一组终端表达式${t_i}$，并使用了一组算符${f_i}$来计算的。 如果$z$的每个子表达式都用一个整数来标记，则可以编写一个抽象程序来计算$z$，方法是计算 一个子表达式序列$\langle y_i\rangle^n_{i = 1}$，并设置$z = y_n$。算法5显示了这样一个程序。

![[0250.jpg]]

每个终端表达式$t_i$是程序的字面常量或输入参数。 这包括了系数，基函数和几何量。 一个非终端子表达式$y_i$是将算符$f_i$应用于先前计算过的表达式序列$\langle y_j\rangle_{j\in \mathcal{J}_i}$的结果，其中$\mathcal{J}_i$是表达式有序序列的记号。注意，根据已计算的子表达式，能产生的相同的$z$值的顺序并不是唯一的。 为了正确起见，我们只要求$j\lt i\\ \forall j\in\mathcal{J}_i$，这样子表达式$y_i$的所有依赖关系都在$y_i$之前计算出来了。 特别是，在此抽象算法中，出于符号上的方便，所有终端表达式都最先编号。

计算z的程序可以表示成图，其中每个表达式$y_i$对应于图顶点。 如果$j \in \mathcal{J}_i$，而$y_i$依赖于$y_j$的值，那么这是有向图的一个从$y_i$到$y_j$的边$e = (i, j)$。 更正式地讲，表示计算$z$的图$G$由一组顶点$V$和一组边$E$组成，它们由以下各项定义：

$$ G = (V, E) \tag{17.65} $$

$$ V = \langle v_i \rangle^n_{i=1} = \langle y_i \rangle^n_{i=1} \tag{17.66} $$

$$ E = {e_k} =\bigcup^n_{i=1} {(i, j)\\ \forall j \in \mathcal{J}_i} \tag{17.67} $$

此图明显是有向的，因为它的依赖关系是有向的。 它是无环的，因为一个表达式只能根据现有的表达式构造。 因此，UFL表达式可以用有向无环图（DAG）表示。 在UFL中，此DAG可用两种方式表示。 在定义表达式时，会建立一个称为表达式树的链表示。 从技术上讲，这仍然是DAG，因为可以在多个子表达式中重用顶点，但是这个表示强调了DAG的树状结构。 另一个表示称为计算图，它与上面$G$的定义非常相似。 此表示对形式编译器最有用。 这两个DAG表示的细节将在下面说明。 它们都将图形中顶点的表示作为表达式对象共享，下面将对其进行说明。

17·6·2 表达式对象

回顾一下算法5的非终端表达式$y_i = f_i(\langle y_j \rangle_{j\in \mathcal{J}_i} )$。 算符$f_i$由表达式对象的类来表示，而表达式$y_i$由此类的实例表示。 在UFL实现中，每个表达式对象都是Expr的某个子类的实例。 Expr类是层次结构的超类，其中包含UFL所支持的所有终端表达式类型和算符类型。 Expr有两个直接的子类，Terminal和Operator，它们将表达式类型层次分为两部分，如图17.2所示。

![[0251.jpg]]

所有表达式对象都被认为是不可变的。 一旦构造了表达式对象，将永远无法对其进行修改。 操作表达式应始终会导致创建新对象。 不可变属性确保表达式对象可以在表达式之间重用和共享，而不会在程序的其他部分产生副作用。 这既减少了内存使用，避免了不必要的对象复制，又简化了对常见子表达式的识别。

在表示$y_i$的Expr对象e上调用e.operands()，会返回一个具有表示$\langle y_j \rangle_{j\in \mathcal{J}_i}$的表达式对象的元组。 请注意，也可应用到没有传出边的终端表达式，并且t.operands()返回空元组。 除了修改表达式对象的操作数外，还可以调用e.reconstruct(operands)用修改后的操作数构造相同类型的新表达式对象，其中operands是表达式对象元组。 如果操作数相同，则此函数返回原始对象，从而允许许多算法节省内存而不会带来其他复杂性。 不变式e.reconstruct(e.operands()) == e应该始终成立。

17·6·3 表达式属性

在第17.4.2节中，讨论了UFL的张量代数和索引记号的功能。 表达式可以是标量或具有任意阶和形状的张量值。 因此，每个表达式对象e的值形状为e.shape()，它是在每个张量轴上维数的整数元组。 标量表达式也有shape ()。 另一个重要的属性是表达式中的一组自由索引，可以使用e.free_indices()作为元组获得。 尽管自由索引没有顺序，但为简单起见，它们以Index实例的元组表示。 因此，元组中的排序没有任何意义。

UFL表达式在引用上是透明的，但有一些例外。 引用透明性是指子表达式可以用其值的另一种表示代替，而无需更改表达式的含义。 这里的关键是，在这种情况下，表达式的值包括张量形状和一组自由索引。 另一个重要点是，函数$f (v)$在某点导数$f'(v)|_{v=g}$，依赖于$v = g$附近的函数值。 这种依赖性的结果是，算符类型在微分时很重要，而不仅是微分变量的当前值。 特别是，Variable不能用其表示的表达式替换，因为diff依赖于Variable实例，而不是有值的表达式。 同样，用某个值替换Coefficient会更改包含关于函数系数求导表达式的含义。

以下示例说明了Variable和diff的这个问题。

# UFL code
e = 0
v = variable(e)
f = sin(v)
g = diff(f, v)

这里的$v$是一个取值为0的变量，但是$\sin(v)$不能简化为0，否则$f$的导数将变成0。 此处的正确结果是$g = \cos(v)$。 打印f和g得到字符串sin(var1(0))和d/d[var1(0)] (sin(var1(0)))。 尝试仅设置$v = e$，看看f和g如何变为零。

17·6·4 树表示

表达式树没有单独的数据结构。 它只是查看表达式结构的一种方式。 任何表达式对象e都可以看作是树的根，其中e.operands()返回其子级。 如果某些子项相等，则它们的出现次数将与表达式中出现的次数相同。 因此很容易遍历树节点（也就是DAG中的vi），但是最终无法直接重用子表达式。 DAG中的边不会明确显示，并且只能通过递归遍历树并选择唯一对象来获得顶点列表。

刚度项$\mathrm{grad}\\ u : \mathrm{grad}\\ v$的表达式树如图17.3所示。 终端表达式u和v没有子节点，$\mathrm{grad}\\ u$项本身由带有两个节点的树表示。 每次将算符应用于某些表达式时，它将返回引用其操作数的新树根。 请注意，用户将在使用该语言时应用的是grad和inner函数，而该图中的名称Grad，Inner和Argument是UFL中用来表示表达式对象的Expr子类名称。 换句话说，由grad(u) 获得的表达式梯度，是一个由Grad(u)表示的表达式，而inner(a, b)则给出一个表达式表示Inner(a, b)。 语言和表示的这种分离仅仅是UFL实现中的一种设计选择。

![[0252.jpg]]

17·6·5 图表示

当以树的形式查看表达式时，所有唯一顶点和边的列表都不是直接可用的。 更加直接地表示DAG可以简化或优化许多算法。 UFL包括从任何表达式构建基于DAG表示的数组的工具（计算图）。 计算图$G = (V, E)$是基于平面数组的数据结构，直接反映了方程（17.65）-（17.67）中图的定义。 此表示可直接访问子表达式之间的依赖关系，并允许在唯一顶点上轻松进行迭代。 该图可通过以下几行轻松构建：

# Python code
# from ufl.algorithms import Graph
from ufl.formatting.graph import Graph
G = Graph(expression)
V, E = G

一个数组（Python列表）V用于存储DAG中的顶点 $\langle v_i\rangle^n_{i=1}$。 对于每个顶点$v_i$，存储一个表达式节点$y_i$来表示。 因此，每个顶点的表达式树也是直接可用的，因为每个表达式节点都是其表达式树的根。 边则存储在数组$E$中，整数元组$(i,j)$表示从$v_i$到$v_j$的边； 也就是说，$v_j$是$v_i$的操作数。 图中的顶点列表是使用深度优先遍历的后序来构建的，这保证了对顶点进行拓扑排序，使得$j\lt i\quad \forall j\in \mathcal{J}_i$。

让我们看一个计算图的例子。 以下代码定义一个简单表达式，然后打印图的顶点和边缘。

# Python code
from ufl import *
cell = triangle
V = FiniteElement("Lagrange", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
c = Constant(cell)
f = Coefficient(V)
e = c*f**2*u*v

# from ufl.algorithms import Graph, partition
from ufl.formatting.graph import Graph, partition
G = Graph(e)
V, E, = G

print("str(e) = %s\n" % str(e))
print("\n".join("V[%d] = %s" % (i, v) for (i, v) in enumerate(V)), "\n")
print("\n".join("E[%d] = %s" % (i, e) for (i, e) in enumerate(E)), "\n")

程序输出的摘录如下所示：

# Generated code
V[0] = v_{-2}
...
V[7] = v_{-1} * c_0 * w_1 ** 2
V[8] = v_{-2} * v_{-1} * c_0 * w_1 ** 2
...
E[6] = (8, 0)
E[7] = (8, 7)

最后两个边，显示了顶点8与顶点7和0的依存关系， 因为$v_8 = v_0 v_7$。 运行代码以查看该代码的完整输出。 尝试更改表达式，然后如上查看图。

从边E可以有效地计算相关的数组。 特别是顶点vi在两个方向上所依赖的顶点索引是有用的：

$$ \begin{aligned}V_{out} &= \langle \mathcal{J}i\rangle^n{i=1} \\ V_{in} &= \langle {j|i \in \mathcal{J}j}\rangle^n{i=1}\end{aligned} \tag{17.68} $$

对任何表达式，这些数组可以很容易构造：

# Python code
Vin = G.Vin()
Vout = G.Vout()

存在类似的函数，用于获取所有传入和传出边E的索引。 由UFL构建的计算图有一个不错的特性：没有两个顶点将代表完全相同的表达式。 在图的构建期间，将子表达式插入哈希映射（Python字典）中即可实现此目的。 一些表达式类对参数进行唯一排序，例如a*b和b*a将成为图中的相同顶点。

当实现某些算法时，表达式节点中的自由索引会使线性化图的解释复杂化，因为具有自由索引的表达式对象不代表一个值，而是代表一组值，自由索引的每个值排列都对应一个值。 一种解决方案是在构造图之前应用expand_indices，它将用自由索引所对应的显式固定索引等效表达式，来替换所有表达式。 但是请注意，展开无法重新获得自由索引。 有关此转换的更多信息，请参见第17.8.3节。

17·6·6 划分

UFL旨在作为形式编译器的前端。 由于最终目标是根据表达式生成代码，因此为代码生成过程提供了一些实用程序。 原则上，只需在顶点上迭代并分别为每个操作生成代码，就可以从其计算图为表达式生成正确的代码，基本上是算法5的镜像。 但是，一个好的形式编译器应该能够产生更好的代码。 UFL提供了实用程序，用于根据子表达式的依赖性将计算图划分为子图（划分），这允许基于形式编译器的正交，可以轻松地将子表达式放置在正确的循环集中。 函数partition实现此功能。 每个划分都由一个简单的顶点索引数组表示，并且每个划分都标记有一组依赖项。 默认情况下，这组依赖使用字符串x，c和v%d来分别表示对空间坐标，特定数量的胞元和形式参数（非系数）的依赖关系。

以下示例代码对上面构建的图进行划分，并根据其依存关系将顶点分组打印。

# Python code
partitions, keys = partition(G)
for deps in sorted(partitions.keys()):
    P = partitions[deps]
    print("The following depends on", tuple(deps))
    for i in sorted(P):
        print("V[%d] = %s" % (i, V[i]))

该程序的输出文本见后。 注意，字面常量2没有依赖。 始终可以在编译时预先计算此划分中的表达式。 常数c_0依赖于每个胞元的变化数据，由依赖集的c代表，而不依赖于空间坐标，因此可以将其置于正交循环之外。 函数w_1及其依赖的表达式还依赖于x表示的空间坐标，因此需要针对每个正交点进行计算。 仅依赖于测试或试探功能的表达式用v%d标记，其中数字是UFL用于区分参数的内部计数器。 请注意，此处的测试和试探函数被标记所依赖的空间坐标，但是不依赖胞元的数量。 这仅适用于在局部参考单元上定义的有限元，在这种情况下，可以在每个正交点中预先计算基函数。 有限元空间中基本函数在运行时的实际依赖关系对于UFL是未知的，这就是为什么函数partition用可选的多功能参数，以便形式编译器可以提供更准确依赖关系的原因。 有关此详细的实现信息，请参考partition的实现。

# Generated code
The following depends on ()
V[4] = 2
The following depends on ("c",)
V[2] = c_0
The following depends on ("x", "c")
V[3] = w_1
V[5] = w_1 ** 2
V[6] = c_0 * w_1 ** 2
The following depends on ("x", "v-1")
V[1] = v_{-1}
The following depends on ("x", "c", "v-1")
V[7] = v_{-1} * c_0 * w_1 ** 2
The following depends on ("x", "v-2")
V[0] = v_{-2}
The following depends on ("x", "c", "v-2", "v-1")
V[8] = v_{-2} * v_{-1} * c_0 * w_1 ** 2

【[[docs/fem/_index#第十七章 UFL：有限元形式语言|章节目录]]】