DOLFIN：C++/Python有限元库》线性代数【翻译】

2021-01-07

| 科学计算 | | FEniCS , 有限元法 , 偏微分方程 |

10·3 功能

DOLFIN被组织为库（模块）的集合，每个库都覆盖某个功能区域。我们在这里回顾这些领域，并解释最常用的类和函数的用途和用法。回顾是自下而上的；也就是说，我们首先描述DOLFIN的核心低层功能（线性代数和网格），然后向上移动以描述高层功能。有关更多详细信息，请参阅FEniCS项目网页上的DOLFIN程序员参考和Logg and Wells（2010）。

10·3·1 线性代数

DOLFIN提供了一系列线性代数对象和函数，包括向量，稠密和稀疏矩阵，直接和迭代线性求解器以及特征值求解器，并通过一个简单且一致的接口来实现。对于大多数基础功能，DOLFIN依赖于第三方库，例如PETSc和Trilinos。 DOLFIN定义了抽象基类GenericTensor，GenericMatrix和GenericVector，它们在整个库中被了广泛使用。在DOLFIN中提供了针对多个后端这些通用接口的实现，从而实现了针对不同后端的通用接口。用户还可以通过实现通用接口来包装其他线性代数后端。

矩阵和向量。创建矩阵和向量的最简单方法是通过Matrix和Vector类。通常，Matrix和Vector表示分布式线性代数对象，当并行运行时，这些对象可以跨（MPI）进程存储。与有限元库中最常用的用法一致，Vector使用稠密存储，而Matrix使用稀疏存储。可以如下创建向量：

// C++ code
Vector x;

# Python code
x = Vector()

可以通过以下方式创建矩阵：

// C++ code
Matrix A;

# Python code
A = Matrix()

在大多数应用程序中，用户可能需要创建矩阵或向量，但是对线性代数对象的大多数操作（包括调整大小）将在库内部进行，并且用户无需直接对对象进行操作。

以下代码说明了如何创建大小为100的向量：【原文的代码已经过期了，被我注释掉了】

// C++ code
// Vector x(100);
x.init(100);

# Python code
# x = Vector(100)
x.init(100)

许多后端支持分布式线性代数用于并行计算，在这种情况下，向量x的全局大小为100，而DOLFIN将在（近）相等大小的部分中跨进程划分向量。

创建一个给定大小的矩阵（稀疏矩阵）要更复杂，因此通常需要根据稀疏矩阵的结构（稀疏模式）进行初始化（分配数据结构）。必要时，稀疏矩阵的初始化由DOLFIN处理。

尽管DOLFIN支持分布式线性代数对象用于并行计算，但很少有用户在并行数据布局级别接触到细节。库将自动处理跨进程的对象分配。

求解线性系统。解决线性系统 $$Ax = b$$ 的最简单方法是使用

// C++ code
solve(A, x, b);

# Python code
solve(A, x, b)

DOLFIN将使用默认方法来求解方程组。可以给出可选的参数来指定求解线性系统时使用的算法，对于迭代方法，则指定使用的前置条件：

// C++ code
solve(A, x, b, "lu");
solve(A, x, b, "gmres", "ilu");

# Python code
solve(A, x, b, "lu");
solve(A, x, b, "gmres", "ilu")

可用的方法和预处理器取决于DOLFIN配置了哪个线性代数后端。要列出可用的求解器方法和预处理器，可以使用以下命令：

// C++ code
list_lu_solver_methods();
list_krylov_solver_methods();
list_krylov_solver_preconditioners();

# Python code
list_lu_solver_methods()
list_krylov_solver_methods()
list_krylov_solver_preconditioners()

使用函数solve很简单，但是对解决过程的细节几乎没有控制。对于许多应用程序，希望对求解过程进行一定程度的控制，并在整个模拟过程中重用求解器对象。

线性系统 $$Ax = b$$ 可以使用LU分解（一种直接方法）求解，如下所示：

// C++ code
LUSolver solver(A);
solver.solve(x, b);

# Python code
solver = LUSolver(A)
solver.solve(x, b)

或者，可以在构造后设置与线性求解器关联的算符 $$A$$ ：

// C++ code
LUSolver solver;
solver.set_operator(A);
solver.solve(x, b);

# Python code
solver = LUSolver()
solver.set_operator(A)
solver.solve(x, b)

当通过函数接口传递线性求解器并将算符设置在函数内部时，这很有用。

在某些情况下，对于给定的矩阵 $$A$$ 和不同的向量 $$b$$ ，或对于不同的 $$A$$ ，但具有相同的非零结构，系统 $$Ax = b$$ 可以求解多次。如果A的非零结构不变，那么可以通知LU求解器来重复求解来提高效率：

// C++ code
solver.parameters["same_nonzero_pattern"] = true;

# Python code
solver.parameters["same_nonzero_pattern"] = True

在 $$A$$ 不变的情况下，可以通过重新使用 $$A$$ 的LU分解来显着减少后续求解时间。分解重用由参数reuse_factorization控制。

某些后端可能会规定要使用的特定LU求解器。这取决于DOLFIN为后端配置了哪些求解器，以及如何配置第三方线性代数后端。

方程组 $$Ax = b$$ 可以使用预处理Krylov求解器通过以下方式求解：

// C++ code
KrylovSolver solver(A);
solver.solve(x, b);

# Python code
solver = KrylovSolver(A)
solver.solve(x, b)

上面将使用默认的预处理器和求解器以及默认参数。如果构造的KrylovSolver没有矩阵算符A，则可以在构造后设置算符：

// C++ code
KrylovSolver solver;
solver.set_operator(A);
solver.solve(x, b);

# Python code
solver = KrylovSolver()
solver.set_operator(A)
solver.solve(x, b)

在某些情况下，使用不同于 $$A$$ 的前置条件矩阵 $$P$$ 可能会有用：

// C++ code
KrylovSolver solver;
solver.set_operators(A, P);
solver.solve(x, b);

# Python code
solver = KrylovSolver()
solver.set_operators(A, P)
solver.solve(x, b)

可以设置Krylov求解器的各种参数。一些常见的参数是：

# Python code
solver = KrylovSolver()
solver.parameters["relative_tolerance"] = 1.0e-6
solver.parameters["absolute_tolerance"] = 1.0e-15
solver.parameters["divergence_limit"] = 1.0e4
solver.parameters["maximum_iterations"] = 10000
solver.parameters["error_on_nonconvergence"] = True
solver.parameters["nonzero_initial_guess"] = False

可以类似地从C++设置参数。可以通过参数控制打印KrylovSolver收敛的摘要和收敛历史的详细信息：可以通过参数控制打印KrylovSolver收敛的摘要和收敛历史的详细信息：

// C++ code
KrylovSolver solver;
solver.parameters["report"] = true;
solver.parameters["monitor_convergence"] = true;

# Python code
solver = KrylovSolver()
solver.parameters["report"] = True
solver.parameters["monitor_convergence"] = True

可以在构造求解器对象时设置要使用的特定Krylov求解器和预处理器。最简单的方法是通过字符串描述设置Krylov方法和前置条件。例如：

// C++ code
KrylovSolver solver("gmres", "ilu");

# Python code
solver = KrylovSolver("gmres", "ilu")

上面指定了通用最小残差（GMRES）方法作为求解器，以及不完全LU（ILU）预处理。

当配置好后端（如PETSc和Trilinos），可以应用各种各样的Krylov方法和预处理器，并且有大量的求解器和预处理器参数可设置。除了此处描述的内容外，DOLFIN还提供了更高级的接口，允许对求解过程进行更精细的控制。用户也可以提供自己的预处理器。

求解特征值问题。 DOLFIN使用基于PETSc的库SLEPc来求解特征值问题。 SLEPc接口仅适用于基于PETSc的线性代数对象。因此，有必要使用基于PETSc的对象，或将默认的线性代数后端设置为PETSc和向下转换的对象（如下一节所述）。以下代码说明了特征值问题 $Ax =\lambda x$ 的解：

// C++ code

// Create matrix
PETScMatrix A;

// Code omitted for setting the entries of A

// Create eigensolver
SLEPcEigenSolver eigensolver(A);

// Compute all eigenvalues of A
eigensolver.solve();

// Get first eigenpair
double lambda_real, lambda_complex;
PETScVector x_real, x_complex;
eigensolver.get_eigenpair(lambda_real, lambda_complex, x_real, x_complex, 0);

# Python code

# Create matrix
A = PETScMatrix()

# Code omitted for setting the entries of A

# Create eigensolver
eigensolver = SLEPcEigenSolver(A)

# Compute all eigenvalues of A
eigensolver.solve()

# Get first eigenpair
lambda_r, lambda_c, x_real, x_complex = eigensolver.get_eigenpair(0)

特征值的实分量和复分量分别以lambda_real和lambda_complex返回，特征向量的实分量和复分量分别以x_real和x_complex返回。

要为广义特征值问题 $Ax =\lambda Mx$ 创建求解器，可以使用 $$A$$ 和 $$M$$ 来构造特征求解器：

// C++ code
PETScMatrix A;
PETScMatrix M;

// Code omitted for setting the entries of A and M
SLEPcEigenSolver eigensolver(A, M);

# Python code
A = PETScMatrix()
M = PETScMatrix()

# Code omitted for setting the entries of A and M

eigensolver = SLEPcEigenSolver(A, M)

用户可以通过参数系统设置许多选项，以控制特征值问题求解过程。要打印可用参数的列表，请分别从C ++和Python调用info(eigensolver.parameters, true)和 info(eigensolver.parameters, True)。

选择线性代数后端。 Matrix，Vector，LUSolver和KrylovSolver对象均基于特定的线性代数后端。默认后端取决于配置DOLFIN时启用了哪些后端。可以通过全局参数linear_algebra_backend设置后端。将PETSc用作线性代数后端：

// C++ code
parameters["linear_algebra_backend"] = "PETSc";

# Python code
parameters["linear_algebra_backend"] = "PETSc"

在创建线性代数对象之前，应先设置此参数。要使用Trilinos集合中的Epetra，应将参数linear_algebra_backend设置为Epetra。对于uBLAS，参数应设置为uBLAS，对于MTL4，参数应设置为MTL4。

用户可以显式创建使用特定后端的线性代数对象。通常，此类对象以后端名称开头。例如，通过以下方法创建基于PETSc的向量和LU求解器：

// C++ code
PETScVector x;
PETScLUSolver solver;

# Python code
x = PETScVector()
solver = PETScLUSolver()

求解非线性系统。 DOLFIN提供了类NewtonSolver的牛顿求解器，用于求解如下的非线性方程组

F(x) = 0 \tag{10.1}

其中， $x \in \mathbb{R}^n$ 和 $F : \mathbb{R}^n \to \mathbb{R}^n$ 。为了使用DOLFIN的牛顿求解器解决此类问题，用户需要提供NonlinearProblem的子类。 NonlinearProblem对象的目的是求 $$F$$ 及其雅可比行列式 $J : \mathbb{R}^n \to \mathbb{R}^n \times \mathbb{R}^n$ 。下面是一个用户提供的用于求解非线性微分方程的MyNonlinearProblem类的概要。

// C++ code
class MyNonlinearProblem : public NonlinearProblem
{
public:

// Constructor
MyNonlinearProblem(const Form& L, const Form& a,
                   const BoundaryCondition& bc) : L(L), a(a), bc(bc) {}

// User-defined residual vector F
void F(GenericVector& b, const GenericVector& x)
{
    assemble(b, L);
    bc.apply(b, x);
}

// User-defined Jacobian matrix J
void J(GenericMatrix& A, const GenericVector& x)
{
    assemble(A, a);
    bc.apply(A);
}

private:
    const Form& L;
    const Form& a;
    const BoundaryCondition& bc;
};

MyNonlinearProblem对象是使用线性形式L（组装时对应于F）和双线性形式a（组装时对应于J）构造的。稍后将在例子中详细讨论使Form和BoundaryCondition类的使用。可以在Python中定义相同的MyNonlinearProblem类：

# Python code
class MyNonlinearProblem(NonlinearProblem):
    def __init__(self, L, a, bc):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bc = bc
    def F(self, b, x):
        assemble(self.L, tensor=b)
        self.bc.apply(b, x)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        self.bc.apply(A)

一旦定义了非线性问题的类，就可以创建一个NewtonSolver对象，并且可以使用牛顿求解器来计算非线性问题的解向量 $$x$$ ：

// C++ code
MyNonlinearProblem problem(L, a, bc);
NewtonSolver newton_solver;

newton_solver.solve(problem, u.vector());

# Python code
problem = MyNonlinearProblem(L, a, bc)
newton_solver = NewtonSolver()

newton_solver.solve(problem, u.vector())

NewtonSolver有许多参数可设。确定牛顿求解器行为的一些参数有：

# Python code
newton_solver = NewtonSolver()
newton_solver.parameters["maximum_iterations"] = 20
newton_solver.parameters["relative_tolerance"] = 1.0e-6
newton_solver.parameters["absolute_tolerance"] = 1.0e-10
newton_solver.parameters["error_on_nonconvergence"] = False

类似地，C++也可设置参数。在测试收敛性时，通常会检查残差F的范数。有时，用检查迭代校正 $$dx$$ 的范数替代会很有用的。这由参数convergence_criterion控制，可以将其设置为"residual"以检查残差F的大小，或设置为"incremental"以检查增量 $$dx$$ 的大小。

对于更高级的用法，可以使用指定在求解过程中使用的线性求解器和前置条件的构造的NewtonSolver。

【章节目录】