3. 热传导方程 The Heat Equation#
本章内容
本章以热传导方程为例介绍时间相关问题的基本程序结构: 连续与离散变分形式、后向 Euler 时间离散、时间循环、收敛阶、时变系数和结果输出. 最后用一个独立示例说明如何移动网格. 后续的 Cahn–Hilliard 和 Navier–Stokes 方程沿用相同的时间循环结构.
Overview
Using the heat equation, this chapter introduces the basic program structure for time-dependent problems: continuous and discrete variational forms, backward Euler time discretization, time loops, convergence rates, time-dependent coefficients, and output. A separate example then illustrates mesh motion. The Cahn–Hilliard and Navier–Stokes chapters use the same time-loop structure.
3.1. 热传导方程与变分形式 The heat equation and its variational forms#
在 \(\Omega=(0,1)^2\) 上考虑
On \(\Omega=(0,1)^2\), consider
其中常数 \(a>0\) 是热扩散率. 记 \((w,v)=\int_\Omega wv\,\mathrm dx\) 为 \(L^2(\Omega)\) 内积, 并定义连续试探和检验空间
Here the constant \(a>0\) is the thermal diffusivity. Let \((w,v)=\int_\Omega wv\,\mathrm dx\) denote the \(L^2(\Omega)\) inner product and define the continuous trial and test spaces
对扩散项分部积分, 得到弱形式: 对几乎处处的 \(t\in(0,T]\), 求 \(u(t)\in V_g(t)\) 使得
Integrating the diffusion term by parts gives the weak form: for almost every \(t\in(0,T]\), find \(u(t)\in V_g(t)\) such that
设 \(\mathcal T_h\) 是 \(\Omega\) 的三角剖分, \(h\) 表示网格尺寸. 定义 \(k\) 次连续 Lagrange 有限元空间及其齐次子空间为
Let \(\mathcal T_h\) be a triangulation of \(\Omega\) and let \(h\) denote its mesh size. The continuous Lagrange finite element space of degree \(k\) and its homogeneous subspace are
非齐次边界值用其离散表示 \(g_h\) 施加, 相应仿射空间记为 \(V_{h,g}(t)\). 将 \([0,T]\) 等分为 \(N\) 个时间步, 记 \(\tau=T/N\) 和 \(t_n=n\tau\). 后向 Euler 格式为: 求 \(u_h^{n+1}\in V_{h,g}(t_{n+1})\), 使得
Nonhomogeneous boundary data are imposed through a discrete representation \(g_h\), with the corresponding affine space denoted by \(V_{h,g}(t)\). Divide \([0,T]\) into \(N\) steps and let \(\tau=T/N\) and \(t_n=n\tau\). At every step, the backward Euler scheme finds \(u_h^{n+1}\in V_{h,g}(t_{n+1})\) such that
3.2. Firedrake 实现 Firedrake implementation#
下面的函数 solve_heat 实现上述后向 Euler 格式, 并由主算例和后面的收敛阶测试共同调用. 为简单起见, 取 \(a=1/(2\pi^2)\)、\(f=g=0\) 和 \(u_0=\sin(\pi x)\sin(\pi y)\). 此时真解为 \(u(\boldsymbol x,t)=e^{-t}\sin(\pi x)\sin(\pi y)\).
The function solve_heat below implements the backward Euler scheme and is reused by both the main example and the convergence tests. For simplicity, take \(a=1/(2\pi^2)\), \(f=g=0\), and \(u_0=\sin(\pi x)\sin(\pi y)\). The exact solution is \(u(\boldsymbol x,t)=e^{-t}\sin(\pi x)\sin(\pi y)\).
代码先将离散方程写成残量形式
The code first writes the discrete equation in residual form
对于一般的源项 \(f\), 双线性形式 \(A(u,v)\) 和线性形式 \(L(v)\) 分别为
For a general source \(f\), the bilinear form \(A(u,v)\) and linear form \(L(v)\) are
本例取 \(f=0\). F 按照 \(A(u_h^{n+1},v_h)-L(v_h)\) 构造, 再由 lhs(F) 和 rhs(F) 自动分离双线性部分和线性部分.
Here \(f=0\). The residual F is formed as \(A(u_h^{n+1},v_h)-L(v_h)\), after which lhs(F) and rhs(F) extract the bilinear and linear parts automatically.
from firedrake import *
from firedrake.pyplot import tricontourf
import matplotlib.pyplot as plt
import numpy as np
def solve_heat(N, num_steps, degree=1, T=1.0):
"""Solve the manufactured heat problem and return its final state."""
dt_value = T/num_steps
dt = Constant(dt_value)
t = Constant(0.0)
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, 'CG', degree)
x, y = SpatialCoordinate(mesh)
a = Constant(1/(2*pi**2))
u_exact = exp(-t)*sin(pi*x)*sin(pi*y)
u_trial = TrialFunction(V)
v_test = TestFunction(V)
u_n = Function(V, name='u_n')
u_h = Function(V, name='u_h')
u_n.interpolate(u_exact)
u_initial = Function(V).assign(u_n)
F = (
(u_trial-u_n)/dt*v_test*dx
+ a*inner(grad(u_trial), grad(v_test))*dx
)
bc = DirichletBC(V, 0, 'on_boundary')
problem = LinearVariationalProblem(lhs(F), rhs(F), u_h, bcs=bc)
solver = LinearVariationalSolver(problem)
for step in range(num_steps):
t.assign((step+1)*dt_value)
solver.solve()
u_n.assign(u_h)
return mesh, u_h, u_exact, u_initial, t
变量
u_trial、v_test: 分别表示新时间层的试探函数 \(u_h^{n+1}\) 和检验函数 \(v_h\).u_n、u_h:u_n保存已知的上一时间层解,u_h接收当前时间步求得的新解.dt、t:dt是时间步长;t是可更新的 UFLConstant, 用于定义随时间变化的u_exact.F: 按“左端减右端”构造的残量, 求解目标是使F = 0.bc、problem、solver: 依次保存边界条件、线性变分问题和可重复使用的求解器.
函数 (算子)
TrialFunction(V)、TestFunction(V): 在离散空间V上创建试探函数和检验函数.lhs(F)、rhs(F): 从线性残量中分别提取双线性形式 \(A\) 和线性形式 \(L\).LinearVariationalProblem(...): 表示求u_h, 使 \(A(u_h,v_h)=L(v_h)\), 同时施加bc.LinearVariationalSolver(problem): 构造线性求解器; 它位于时间循环外, 因而可以在各时间步重复使用.t.assign(...): 更新时间系数;u_exact会读取新值, 不必重建表达式.solver.solve()与u_n.assign(u_h): 前者计算新解, 后者随即把新解保存为下一步的旧解.solve_heat(...): 封装完整计算并返回验证所需对象, 后面的收敛阶测试可以直接调用.
对非线性问题通常直接把残量交给 NonlinearVariationalProblem, 不使用 lhs 和 rhs.
Variables
u_trial,v_test: the trial function \(u_h^{n+1}\) at the new time level and the test function \(v_h\), respectively.u_n,u_h:u_nstores the known solution from the previous time level, whileu_hreceives the newly computed solution.dt,t:dtis the time-step size;tis an updatable UFLConstantused in the time-dependent expressionu_exact.F: the residual formed as “left-hand side minus right-hand side”; the goal is to enforceF = 0.bc,problem,solver: store the boundary condition, linear variational problem, and reusable solver, respectively.
Functions (operators)
TrialFunction(V),TestFunction(V): create the trial and test functions on the discrete spaceV.lhs(F),rhs(F): extract the bilinear form \(A\) and linear form \(L\) from the linear residual.LinearVariationalProblem(...): asks foru_hsatisfying \(A(u_h,v_h)=L(v_h)\) while applyingbc.LinearVariationalSolver(problem): constructs the linear solver outside the time loop so that it can be reused at every step.t.assign(...): updates the time coefficient;u_exactreads the new value without rebuilding the expression.solver.solve()andu_n.assign(u_h): the former computes the new solution, and the latter stores it as the old solution for the next step.solve_heat(...): encapsulates the complete computation and returns the objects needed for verification, so the convergence tests can call it directly.
Nonlinear problems normally pass the residual directly to NonlinearVariationalProblem instead of using lhs and rhs.
下面在固定网格上求解该问题, 并计算 \(L^2\) 误差.
We now solve the problem on a fixed mesh and compute the \(L^2\) error.
mesh, u_h, u_exact, u_initial, t = solve_heat(
N=32, num_steps=32, degree=1, T=1.0
)
error_L2 = sqrt(assemble((u_exact-u_h)**2*dx(domain=mesh, degree=6)))
print(f'final time = {float(t):.2f}, L2 error = {error_L2:.3e}')
final time = 1.00, L2 error = 2.118e-03
上面输出终止时刻的 \(L^2\) 误差, 用于定量比较离散解与解析解; 其收敛行为将在下一节通过加密时间步和网格来检验. 下图使用相同的色标比较初值与终值; 解析解的空间形状不变, 振幅按 \(e^{-t}\) 衰减.
The final-time \(L^2\) error is printed for a quantitative comparison between the discrete and exact solutions. Its convergence behavior is examined in the next section by refining the time step and the mesh. The plots use a common color scale to compare the initial and final states; the exact solution keeps the same spatial shape while its amplitude decays as \(e^{-t}\).
fig, axes = plt.subplots(1, 2, figsize=(8, 3.3), constrained_layout=True)
for ax, uh, title in zip(axes, (u_initial, u_h), ('$t=0$', f'$t={float(t):.0f}$')):
colors = tricontourf(uh, axes=ax, levels=np.linspace(0, 1, 11))
ax.set_aspect('equal')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_title(title)
_ = fig.colorbar(colors, ax=ax, shrink=0.82)
3.3. 时间与空间收敛阶 Temporal and spatial convergence rates#
对光滑解, 上述数值格式有如下误差估计
For a smooth solution, the scheme above satisfies the error estimates
验证时间阶时固定网格并使用四次元, 然后逐次减半 \(\tau\). 验证空间阶时使用一次元并令 \(\tau=h^2\). 这里取较短的终止时间 \(T=1/4\) 以控制测试开销. solve_and_measure 只调用前面的 solve_heat, 再用 compute_errors 计算误差.
For the temporal test, we fix the mesh, use degree-four elements, and halve \(\tau\) successively. For the spatial test, we use linear elements and set \(\tau=h^2\). The shorter final time \(T=1/4\) keeps the cost of the tests under control. solve_and_measure simply calls solve_heat above and evaluates the errors with compute_errors.
def compute_errors(mesh, u_h, u_exact):
error = u_exact-u_h
dx_error = dx(domain=mesh, degree=10)
error_L2 = sqrt(assemble(error**2*dx_error))
error_H1 = sqrt(assemble((error**2+inner(grad(error), grad(error)))*dx_error))
return error_L2, error_H1
def solve_and_measure(N, num_steps, degree):
mesh, u_h, u_exact, _, _ = solve_heat(
N, num_steps, degree=degree, T=0.25
)
return compute_errors(mesh, u_h, u_exact)
def convergence_rates(errors, scales):
rates = [np.nan]
for i in range(1, len(errors)):
rates.append(np.log(errors[i-1]/errors[i])/np.log(scales[i-1]/scales[i]))
return rates
def print_convergence_table(title, headers, rows):
rows = [[str(value) for value in row] for row in rows]
widths = [
max(len(header), *(len(row[j]) for row in rows))
for j, header in enumerate(headers)
]
print(title)
print(' '.join(header.rjust(width) for header, width in zip(headers, widths)))
for row in rows:
print(' '.join(value.rjust(width) for value, width in zip(row, widths)))
3.3.1. 时间收敛阶 Temporal convergence rate#
时间步数依次加倍, 即 \(\tau\) 依次减半. 空间采用固定的 \(8\times8\) 网格和四次元.
The number of time steps is doubled successively, so \(\tau\) is halved each time. Space is discretized on a fixed \(8\times8\) mesh with degree-four elements.
time_steps = [4, 8, 16, 32]
time_dts = [0.25/M for M in time_steps]
time_errors = [solve_and_measure(8, M, 4)[0] for M in time_steps]
time_rates = convergence_rates(time_errors, time_dts)
time_rows = [
[steps, f'{dt_i:.3e}', f'{error_i:.3e}',
'--' if np.isnan(rate_i) else f'{rate_i:.2f}']
for steps, dt_i, error_i, rate_i
in zip(time_steps, time_dts, time_errors, time_rates)
]
print_convergence_table(
'Temporal convergence', ['steps', 'dt', 'L2 error', 'rate'], time_rows
)
assert time_rates[-1] > 0.9
Temporal convergence
steps dt L2 error rate
4 6.250e-02 2.932e-03 --
8 3.125e-02 1.493e-03 0.97
16 1.562e-02 7.535e-04 0.99
32 7.812e-03 3.785e-04 0.99
3.3.2. 空间收敛阶 Spatial convergence rates#
对一次元逐次加密网格, 并同步取 \(\tau=h^2\). 表中分别给出 \(L^2\) 和 \(H^1\) 误差及对应收敛阶.
We refine the mesh successively for linear elements and simultaneously set \(\tau=h^2\). The table reports both \(L^2\) and \(H^1\) errors and their corresponding rates.
如何配合时间步长与网格尺寸 Balancing the time step and mesh size
更一般地, 假设误差估计为
More generally, suppose the error estimate is
并且已经验证时间收敛阶为 \(q\). 测试空间收敛阶时可令 \(\tau=h^s\), 此时
and the order \(q\) with respect to \(\tau\) has already been verified. To test the order with respect to \(h\), set \(\tau=h^s\). Then
应取 \(sq\geq p\), 使时间误差不会降低关于 \(h\) 的观测阶; 若希望时间误差比空间误差渐近地更小, 则取 \(sq>p\).
Choose \(sq\geq p\) so that the temporal error does not reduce the observed order in \(h\); if the temporal error should be asymptotically smaller than the spatial error, choose \(sq>p\).
space_Ns = [8, 16, 32]
space_hs = [1/N for N in space_Ns]
space_steps = [round(0.25/h**2) for h in space_hs]
space_errors = [solve_and_measure(N, M, 1) for N, M in zip(space_Ns, space_steps)]
space_L2 = [errors[0] for errors in space_errors]
space_H1 = [errors[1] for errors in space_errors]
space_rates_L2 = convergence_rates(space_L2, space_hs)
space_rates_H1 = convergence_rates(space_H1, space_hs)
space_rows = [
[N, steps, f'{h_i:.3e}', f'{0.25/steps:.3e}',
f'{e0:.3e}', '--' if np.isnan(r0) else f'{r0:.2f}',
f'{e1:.3e}', '--' if np.isnan(r1) else f'{r1:.2f}']
for N, steps, h_i, e0, r0, e1, r1 in zip(
space_Ns, space_steps, space_hs, space_L2,
space_rates_L2, space_H1, space_rates_H1
)
]
print_convergence_table(
'Spatial convergence',
['N', 'steps', 'h', 'dt', 'L2 error', 'rate', 'H1 error', 'rate'],
space_rows,
)
assert space_rates_L2[-1] > 1.8
assert space_rates_H1[-1] > 0.9
Spatial convergence
N steps h dt L2 error rate H1 error rate
8 16 1.250e-01 1.562e-02 1.491e-02 -- 3.368e-01 --
16 64 6.250e-02 3.906e-03 3.790e-03 1.98 1.695e-01 0.99
32 256 3.125e-02 9.766e-04 9.516e-04 1.99 8.488e-02 1.00
3.4. 时间循环、时变系数与结果输出 Time loops, time-dependent coefficients, and output#
时间相关问题通常都采用同一个程序框架: 在循环外构造变分问题、边界条件和求解器; 在循环内更新时间或其他系数, 求新时间层解, 再更新保存的旧解. 上面的 solve_heat 已经展示了线性问题的基本写法. 线性问题可以直接复用求解器; 对以后将介绍的 Cahn–Hilliard 和 Navier–Stokes 等非线性问题, 还常用旧解初始化当前未知函数, 以减少 Newton 迭代次数.
若源项、边界值或材料参数随时间变化, 应将时间定义为 Constant, 再用包含它的 UFL 表达式定义系数. 时间循环中调用 t.assign(...) 后, 已构造的变分形式会读取新值, 不必重新构造求解器或重新编译内核. 上面的真解 u_exact 就采用了这种写法.
需要查看完整时间序列时, 可以在循环前创建 VTKFile, 然后按适当间隔写出结果. time= 参数使 ParaView 能按物理时间播放数据; 不必每一步都写出. 并行运行时仍使用同一写法, Firedrake 会负责生成分区数据文件.
Time-dependent problems normally share the same program structure: construct the variational problem, boundary conditions, and solver outside the loop; update time or other coefficients inside the loop, solve for the new time level, and then update the stored old solution. The solve_heat function above has already illustrated the basic pattern for a linear problem, whose solver can be reused directly. For nonlinear problems such as the Cahn–Hilliard and Navier–Stokes equations introduced in later chapters, the old solution is also commonly used to initialize the current unknown and reduce the number of Newton iterations.
If a source term, boundary value, or material parameter varies in time, define time as a Constant and use it in the corresponding UFL expression. After t.assign(...) in the loop, the already constructed form reads the new value; neither the solver nor the compiled kernel needs to be rebuilt. The exact solution u_exact above follows this pattern.
To inspect the complete time series, create a VTKFile before the loop and write the solution at a suitable interval. The time= argument lets ParaView play the data in physical time, and output need not be written at every step. The same code works in parallel, with Firedrake generating the partitioned data files.
output = VTKFile('result.pvd')
output.write(u_initial, time=0.0)
output_interval = 10
for step in range(num_steps):
# 更新时间、求解并更新旧解
...
if (step+1) % output_interval == 0:
output.write(u_h, time=float(t))
在 ParaView 中, Warp By Scalar 可显示标量曲面, Glyph 可显示向量箭头, Stream Tracer 可生成流线, Connectivity 可区分并行计算中的不同连通区域或分区. 对当前书中的二维标量结果, 通常先在网页中用 tricontourf; 只有需要交互式观察时间变化或三维数据时才写出 PVD 文件.
In ParaView, Warp By Scalar displays a scalar surface, Glyph displays vector arrows, Stream Tracer generates streamlines, and Connectivity can distinguish connected regions or partitions in parallel data. For the two-dimensional scalar results in this book, tricontourf in the web page is usually sufficient; PVD output is most useful for interactive time-dependent or three-dimensional visualization.
3.5. 移动网格 Moving meshes#
3.5.1. 物质导数、ALE 弱形式与离散格式 Material derivative, ALE weak form, and discretization#
一个简单的移动边界问题是在膨胀圆盘 \(\Omega(t)=\{\boldsymbol x:|\boldsymbol x|<e^t\}\) 上求解热传导方程
A simple moving-boundary problem is to solve the heat equation on the expanding disk \(\Omega(t)=\{\boldsymbol x:|\boldsymbol x|<e^t\}\):
设 ALE 映射 \(\boldsymbol\chi(\boldsymbol X,t):\Omega(0)\to\Omega(t)\) 描述网格点的运动, 并定义网格速度
沿网格轨迹观察 \(u\) 得到相对于网格速度的物质导数 (也称 ALE 导数)
所以 \(\partial_tu=D_t^{\boldsymbol w}u-\boldsymbol w\cdot\nabla u\). 记 \(V_g(t)=\{v\in H^1(\Omega(t)):v|_{\partial\Omega(t)}=g\}\) 和 \(V_0(t)=H_0^1(\Omega(t))\). 在当前区域上, ALE 弱形式为: 求 \(u(t)\in V_g(t)\), 使对任意随网格输运的 \(v(t)\in V_0(t)\),
Let the ALE map \(\boldsymbol\chi(\boldsymbol X,t):\Omega(0)\to\Omega(t)\) describe the motion of the mesh points, and define the mesh velocity
Observing \(u\) along the mesh trajectories gives the material derivative relative to the mesh velocity (also called the ALE derivative)
so \(\partial_tu=D_t^{\boldsymbol w}u-\boldsymbol w\cdot\nabla u\). Let \(V_g(t)=\{v\in H^1(\Omega(t)):v|_{\partial\Omega(t)}=g\}\) and \(V_0(t)=H_0^1(\Omega(t))\). On the current domain, the ALE weak form reads: find \(u(t)\in V_g(t)\) such that, for every mesh-transported \(v(t)\in V_0(t)\),
令 \(\Phi^n:\Omega^n\to\Omega^{n+1}\) 为一步网格映射, \(\widetilde u_h^n=u_h^n\circ(\Phi^n)^{-1}\) 为把旧解随网格搬到新区域后的函数, \(\boldsymbol w_h^{n+1}\) 为该步的离散网格速度. 后向 Euler 格式为: 求 \(u_h^{n+1}\in V_g^{n+1}\), 使对任意 \(v_h\in V_0^{n+1}\),
Let \(\Phi^n:\Omega^n\to\Omega^{n+1}\) be the one-step mesh map, let \(\widetilde u_h^n=u_h^n\circ(\Phi^n)^{-1}\) be the old solution carried to the new domain along with the mesh, and let \(\boldsymbol w_h^{n+1}\) be the discrete mesh velocity of this step. The backward Euler scheme reads: find \(u_h^{n+1}\in V_g^{n+1}\) such that, for every \(v_h\in V_0^{n+1}\),
当 Firedrake 网格坐标原位更新时, Function 的自由度值不会随坐标改变; 因此更新坐标后保留的 u_n 正好表示 \(\widetilde u_h^n\).
When Firedrake mesh coordinates are updated in place, the degree-of-freedom values of a Function do not change with the coordinates; hence the u_n retained after the coordinate update is exactly \(\widetilde u_h^n\).
3.5.2. 最小 ALE 示例 A minimal ALE example#
取 \(\boldsymbol\chi(\boldsymbol X,t)=e^t\boldsymbol X\), 则 \(\Omega(t)\) 是半径为 \(e^t\) 的圆盘, \(\boldsymbol w=\boldsymbol x\). 选择解析解 \(u=e^{2t}-|\boldsymbol x|^2\) 和扩散率 \(a=0.1\), 则零 Dirichlet 边界条件成立, 源项为 \(f=2e^{2t}+4a\). 每一步先按精确 ALE 映射放大网格, 再由坐标增量计算该步的离散网格速度并求解上面的离散方程.
Take \(\boldsymbol\chi(\boldsymbol X,t)=e^t\boldsymbol X\), so \(\Omega(t)\) is the disk of radius \(e^t\) and \(\boldsymbol w=\boldsymbol x\). With diffusivity \(a=0.1\), choose the exact solution \(u=e^{2t}-|\boldsymbol x|^2\), which satisfies the homogeneous Dirichlet condition with source \(f=2e^{2t}+4a\). At each step, the mesh is first expanded by the exact ALE map; the discrete mesh velocity is then computed from the coordinate increment before solving the discrete equation above.
mesh_ale = UnitDiskMesh(3)
V_ale = FunctionSpace(mesh_ale, 'CG', 1)
x_ale = SpatialCoordinate(mesh_ale)
num_ale_steps = 10
dt_ale_value = 0.01
dt_ale = Constant(dt_ale_value)
t_ale = Constant(0.0)
diffusivity = Constant(0.1)
u_ale = TrialFunction(V_ale)
v_ale = TestFunction(V_ale)
u_ale_n = Function(V_ale, name='u')
u_ale_h = Function(V_ale, name='u')
u_ale_exact = exp(2*t_ale)-dot(x_ale, x_ale)
u_ale_n.interpolate(u_ale_exact)
coordinate_space = mesh_ale.coordinates.function_space()
old_coordinates = Function(coordinate_space)
mesh_velocity = Function(coordinate_space)
source = 2*exp(2*t_ale)+4*diffusivity
a_ale = (
u_ale/dt_ale*v_ale
- dot(mesh_velocity, grad(u_ale))*v_ale
+ diffusivity*inner(grad(u_ale), grad(v_ale))
)*dx
L_ale = (u_ale_n/dt_ale+source)*v_ale*dx
bc_ale = DirichletBC(V_ale, 0, 'on_boundary')
problem_ale = LinearVariationalProblem(a_ale, L_ale, u_ale_h, bcs=bc_ale)
solver_ale = LinearVariationalSolver(problem_ale)
for step in range(num_ale_steps):
t_ale.assign((step+1)*dt_ale_value)
old_coordinates.assign(mesh_ale.coordinates)
mesh_ale.coordinates.assign(exp(dt_ale_value)*old_coordinates)
mesh_velocity.assign((mesh_ale.coordinates-old_coordinates)/dt_ale)
solver_ale.solve()
u_ale_n.assign(u_ale_h)
error_ale = sqrt(assemble((u_ale_h-u_ale_exact)**2*dx(degree=6)))
print(f'final time = {float(t_ale):.2f}, ALE L2 error = {error_ale:.3e}')
final time = 0.10, ALE L2 error = 6.686e-03
Firedrake 与 DMPlex 中的坐标
mesh.coordinates.assign(...) 更新的是 Firedrake 坐标函数. 若后续直接调用依赖 DMPlex 坐标的 PETSc 操作, 还需要同步坐标; 具体方法见 更新 DMPlex 中的网格坐标 Updating mesh coordinates in DMPlex.
Coordinates in Firedrake and DMPlex
mesh.coordinates.assign(...) updates Firedrake’s coordinate function. If a later PETSc operation reads coordinates directly from DMPlex, the coordinates must also be synchronized; see 更新 DMPlex 中的网格坐标 Updating mesh coordinates in DMPlex.