复数值有限元方法 — Helmholtz¶
一个脚本 examples/wave/helmholtz/helmholtz.py 在单位正方形上求解带复数值系数的时谐(频域)Helmholtz 方程。它端到端地验证了在 ROADMAP item 2 中解锁的 复数装配路径:复数单元矩阵现在能够从复数 point_data 系数出发,一路流经装配与狄利克雷静态凝聚,进入复数线性求解 —— 并通过正确的复数伴随回传。
强形式是内部 Helmholtz 问题
完全由边界数据驱动。我们采用 人工构造的平面波解
因此狄利克雷数据为 \(g = u_\text{exact}|_{\partial\Omega}\)。体力恰好为零,因为 \(-\Delta e^{ikx} = k^2 e^{ikx}\) 逐点抵消了 \(-k^2 u\) 质量项 —— 这使得解析解在整个域上都可用于误差检验。
为什么复数示例很重要¶
装配栈过去处处假设实数 dtype。Helmholtz 是端到端演练复数路径的最小问题:
复数系数经由
point_data。波数项 \(k^2\) 作为(此处为常数、可为复数的)逐节点场携带,并在ElementAssembler.__call__内部沿单元与求积广播 —— 这正是 PML 层用于其各向异性、随空间变化的复数系数的同一机制。复数狄利克雷静态凝聚。 边界值 \(g\) 是复数,因此
Condenser必须携带复数的内部系统右端项。复数线性求解。
SparseMatrix.solve委托给torch-sla的复对称 \(LDL^\top\) / 厄米 \(LDL^\mathsf{H}\) 分解,它们同时提供匹配的复数伴随 —— 这对下游基于梯度的设计至关重要。
TensorMesh 设置¶
弱形式是单个双线性装配器;唯一标记其为复数的,是 dtype 转换和复数 k_sq 系数:
class HelmholtzAssembler(ElementAssembler):
# a(u, v) = ∫ ∇u·∇v - k² u v dΩ
def forward(self, gradu, gradv, u, v, k_sq):
return gradu @ gradv - k_sq * u * v
mesh = gen_rectangle(chara_length=h, element_type="tri")
points = mesh.points.to(torch.float64).to(device)
# k² as a complex per-node coefficient (constant here; a PML layer
# would make it anisotropic and spatially varying).
k_sq_field = torch.full((mesh.n_points,), k * k + 0j,
dtype=torch.complex128, device=device)
asm = HelmholtzAssembler.from_mesh(mesh, quadrature_order=3)
asm.type(torch.complex128).to(device) # cast the assembler complex
H = asm(points=points, point_data={"k_sq": k_sq_field})
# Dirichlet g = u_exact on the boundary; the inner RHS is complex.
g = torch.exp(1j * k * mesh.points[:, 0].to(torch.float64)).to(torch.complex128)
condenser = Condenser(mesh.boundary_mask, dirichlet_value=g[mesh.boundary_mask])
H_inner, rhs_inner = condenser(H, torch.zeros(mesh.n_points, dtype=torch.complex128))
u = condenser.recover(H_inner.solve(rhs_inner))
两点值得注意:
网格保持实数。 几何、形函数和求积权重保持
float64;只有 系数 和由此得到的系统是复数。asm.type(torch.complex128)将装配器提升为复数,实数形函数张量在装配einsum内部按需上转换到系数的复数 dtype。``complex128`` 是默认值。 收敛性研究在双精度复数下运行;``--dtype complex64`` 可用,但与实数有限元一样,推荐使用双精度以获得干净的收敛率。
收敛性¶
细化网格使 \(L^2\) 误差相对 \(u_\text{exact}\) 以期望的有限元速率下降(在中等 \(k\) 下会有通常的 Helmholtz“污染”)。在 \(k = 2\pi\) 时:
h=0.200 n_dofs= 44 L2 err = 1.529e-01
h=0.100 n_dofs= 143 L2 err = 5.274e-02
h=0.050 n_dofs= 509 L2 err = 1.506e-02
h=0.025 n_dofs=1934 L2 err = 3.935e-03
图 42 helmholtz.py 在 \(k = 2\pi\)、\(h = 0.1\) 时的输出:计算所得场的实部、虚部,以及逐点误差 \(|u - u_\text{exact}|\)。平面波 \(e^{ikx}\) 沿 \(x\) 传播;误差面板在整个域上保持在离散化误差下限。¶
与 scikit-fem 的交叉验证¶
复数路径最强的正确性信号是一条独立的流水线。tests/assemble/test_helmholtz_example.py 把 相同 的 (points, cells) 交给 scikit-fem 的 MeshTri,用其内置的 laplace / mass 积分子装配同一弱形式,以 scipy.sparse.linalg.spsolve 求解,并逐节点比较。在 \(h = 0.1\)、\(k = 2\pi\) 时,两个求解器吻合到浮点精度:
即机器 \(\varepsilon\) —— 两条流水线相对解析平面波继承了 相同 的 \(5.27\times10^{-2}\) 离散化误差。
运行示例¶
cd examples/wave/helmholtz
python helmholtz.py # k = 2π, writes helmholtz.png
python helmholtz.py --k 12.566 --chara-length 0.05 # k = 4π
python helmholtz.py --no-plot # convergence table only