5.1.1. 示例 1: 一维热传导方程

(1) 控制方程

无量纲控制方程为 Eq.5.1.1, \(x \in [0,1]\)。

(5.1.1)\[\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}\]

初始条件为单个方波，即：

(5.1.2)\[\begin{split}u(x, t=0) = \left\{ \begin{array}{lr} 0, & x \le 0.3\\ 1, & x \in (0.3,0.5) \\ 0, & x \ge 0.5\\ \end{array} \right.\end{split}\]

边界条件为零阶条件 (Dirichlet boundary conditions)，即：

(5.1.3)\[u(x=0, t) = u(x=1, t) = 0\]

(2) 离散方法

有限差分方法
均匀结构网格
全离散方法
一阶隐式时间推进
二阶中心差分空间离散

离散方程：

(5.1.4)\[\frac{\hat u_i - u_i}{\Delta t} = \frac{\hat u_{i+1} - \hat u_{i} + \hat u_{i-1}} {\Delta x^2}\]

其中，\(\hat u\) 为下一时间步的解， \(\Delta t, \Delta x\) 为离散单元，\(i\) 为网格点编号。

那么，由 Eq.5.1.4 可得隐式时间推进的单时间步方程为：

(5.1.5)\[\hat u_i + \Delta t / \Delta x^2 (\hat u_{i+1} - \hat u_{i} + \hat u_{i-1}) = u_i\]

则 Eq.5.1.5 可以写为线性方程组形式 \(A\mathbf{u}=\mathbf{b}\), 其中 \(A\) 是一个近似三对角形式的方阵， \(\mathbf{u}\) 是下一时间步的未知解向量， \(\mathbf{b}\) 是当前时间步的解向量。

(3) 数值求解模块

使用 petsc4py 求解线性方程组 Eq.5.1.5, PETSc 的线性求解器可以在运行时指定：

# Use this flag to monitor the solution process:
python example_01.py -ksp_monitor

# Use this flag to swap the linear solver
python example_01.py -ksp_type [SOLVER_NAME]

(4) 代码示例

以下为完整代码。

import petsc4py
import sys
import numpy as np
import matplotlib.pyplot as plt

petsc4py.init(sys.argv)

from petsc4py import PETSc

N_POINTS = 1001
TIME_STEP_LENGTH = 0.005
N_TIME_STEPS = 10

def main():
    
    element_length = 1.0 / (N_POINTS - 1)
    mesh = np.linspace(0.0, 1.0, N_POINTS)

    # Create a new sparse PETSc matrix, fill it and then assemble it
    A = PETSc.Mat().createAIJ([N_POINTS, N_POINTS])
    A.setUp()

    diagonal_entry = 1.0 + 2.0 * TIME_STEP_LENGTH / element_length**2
    off_diagonal_entry = - 1.0 * TIME_STEP_LENGTH / element_length**2

    A.setValue(0, 0, 1.0)
    A.setValue(N_POINTS-1, N_POINTS-1, 1.0)

    for i in range(1, N_POINTS - 1):
        A.setValue(i, i, diagonal_entry)
        A.setValue(i, i-1, off_diagonal_entry)
        A.setValue(i, i+1, off_diagonal_entry)
    
    A.assemble()

    # Define the initial condition
    initial_condition = np.where(
        (mesh > 0.3) & (mesh < 0.5),
        1.0,
        0.0,
    )

    # Assemble the initial rhs to the linear system
    b = PETSc.Vec().createSeq(N_POINTS)
    b.setArray(initial_condition)
    b.setValue(0, 0.0)
    b.setValue(N_POINTS-1, 0.0)

    # Allocate a PETSc vector storing the solution to the linear system
    x = PETSc.Vec().createSeq(N_POINTS)

    # Instantiate a linear solver: Krylow subspace linear iterative solver
    ksp = PETSc.KSP().create()
    ksp.setOperators(A)
    ksp.setFromOptions()

    chosen_solver = ksp.getType()
    print(f"Solving with {chosen_solver:}")
    
    plt.plot(mesh, initial_condition)

    for iter in range(N_TIME_STEPS):
        ksp.solve(b, x)

        # Re-assemble the rhs to move forward in time
        current_solution = x.getArray()
        b.setArray(current_solution)
        b.setValue(0, 0.0)
        b.setValue(N_POINTS - 1, 0.0)

        # Visualize
        plt.plot(mesh, current_solution)        

    plt.savefig('../figures/example_01.jpg', dpi=300)


if __name__ == "__main__":
    main()

(5) 结果展示

../../_images/example_01.jpg — 图 5.1.1 多个时间步的温度分布