复数值Poisson方程#
原文：https://jsdokken.com/dolfinx-tutorial/chapter1/complex_mode.html
作者: Jørgen S. Dokken
翻译：马鹏飞
许多偏微分方程要在复数域上求解，例如Helmholtz 方程。简单起见，考虑如下形式的Poisson方程：
\[-\Delta u = f \text{ in } \Omega,\]
\[ f = -1 - 2j \text{ in } \Omega,\]
\[ u = u_{exact} \text{ on } \partial\Omega,\]
\[u_{exact}(x, y) = \frac{1}{2}x^2 + 1j\cdot y^2,\]
如章节 求解Poisson方程 所述，我们先给出偏微分方程的弱形式。
定义离散的函数空间 \(V_h\)，令\(u_h\in V_h\)。 令\(u_h = \sum_{i=1}^N c_i \phi_i(x, y)\)，其中\(\phi_i\) 是空间\(V_h\)的 实数值 全局基函数，\(c_i \in \mathcal{C}\) 是复数值的自由度。
接下来，我们选择测试函数 \(v\in \hat V_h\)。如前面章节所述， \(\hat V_h\subset V_h\) 且 \(v\vert_{\partial\Omega}=0\)。
我们需要定义我们的内积空间，我们选取的\(L^2\)内积空间定义了*半双线性型*，映射 \(\langle u, v\rangle\)的原象属于\(V_h\times V_h\mapsto K\)，\(\langle u, v \rangle = \int_\Omega u \cdot \bar v ~\mathrm{d} x\)。此半双线性型具有如下性质：
译者的理解：半双线性型为复数值泛函，是双线性型的一般化。
\[\langle u , v \rangle = \overline{\langle v, u \rangle},\]
\[\langle u , u \rangle \geq 0.\]
现在我们在内积空间中进行分部积分。
\[\int_\Omega \nabla u_h \cdot \nabla \overline{v}~ \mathrm{dx} = \int_{\Omega} f \cdot \overline{v} ~\mathrm{d} s \qquad \forall v \in \hat{V}_h.\]
安装支持复数版本的 FEniCSx#
FEniCS支持实数和复数，函数空间中基函数的系数可以是实数值，也可以是复数值。
import dolfinx
from mpi4py import MPI
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
u_r = dolfinx.fem.Function(V, dtype=np.float64)
u_r.interpolate(lambda x: x[0])
u_c = dolfinx.fem.Function(V, dtype=np.complex128)
u_c.interpolate(lambda x:0.5*x[0]**2 + 1j*x[1]**2)
print(u_r.x.array.dtype)
print(u_c.x.array.dtype)
float64
complex128
然而，在求解线性代数问题\(Ax=b\)时，我们需要使用支持实数和复数的矩阵和向量。我们使用的线性代数求解器是 PETSc ，因此程序需要和它对接。
只可惜PETSc的矩阵只支持一种浮点数据类型，因此我们得安装两个版本的PETSc，一个版本支持float64，另一种支持complex128。在DOLFINx的Dokcer镜像里，两种版本都安装了，可以通过命令source dolfinx-real-mode 或 source dolfinx-complex-mode更改。dolfinx/lab镜像可以通过改变Python内核版本来切换，在浏览器的Jupyter页面中，点击 Kernel->Change Kernel... ，通过选择 Python3 (ipykernel) 或 Python3 (DOLFINx complex) 来确定版本。
检查当前的PETSc安装版本。（注：'c'应该是complex，复数的意思。）
from petsc4py import PETSc
print(PETSc.ScalarType)
assert np.dtype(PETSc.ScalarType).kind == 'c'
<class 'numpy.float64'>
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Cell In[2], line 3
1 from petsc4py import PETSc
2 print(PETSc.ScalarType)
----> 3 assert np.dtype(PETSc.ScalarType).kind == 'c'
AssertionError:
变分问题#
接下来，我们可以定义变分问题了。
import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = dolfinx.fem.Constant(mesh, PETSc.ScalarType(-1 - 2j))
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx
我们使用 PETSc.ScalarType 来封装右端常量的源项，这能让我们最后的积分结果是正确的浮点类型
其次，注意到我们使用 ufl.inner 来描述 \(f\) 和 \(v\) 的乘积(当它们是标量时也是成立的)。根据前面定义的 \(L^2\) 内积， ufl.inner 会将第二个参数取共轭。
内积和求导#
L2 = f * ufl.conj(v) * ufl.dx
print(L)
print(L2)
{ c_0 * (conj((v_0))) } * dx(<Mesh #0>[everywhere], {})
{ c_0 * (conj((v_0))) } * dx(<Mesh #0>[everywhere], {})
类似的，如果我们想要使用函数ufl.derivative计算泛函的导数，我们需要特别注意。因为derivative将一个ufl.TestFunction作为变量，我们需要使用它的共轭才能组装向量。这里同样使用了ufl.conj(v)。
J = u_c**2 * ufl.dx
F = ufl.derivative(J, u_c, ufl.conj(v))
residual = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F))
print(residual.array)
[4.21666667e-03+3.33333333e-05j 1.58750000e-03+8.33333333e-06j
4.68333333e-03+1.00000000e-04j 8.13333333e-03+2.66666667e-04j
3.35000000e-03+3.33333333e-05j 4.68333333e-03+3.66666667e-04j
6.43333333e-03+2.66666667e-04j 8.13333333e-03+8.66666667e-04j
2.58333333e-03+3.33333333e-05j 4.68333333e-03+8.33333333e-04j
4.93333333e-03+2.66666667e-04j 6.43333333e-03+8.66666667e-04j
8.13333333e-03+1.86666667e-03j 1.91666667e-03+3.33333333e-05j
4.68333333e-03+1.50000000e-03j 3.63333333e-03+2.66666667e-04j
4.93333333e-03+8.66666667e-04j 6.43333333e-03+1.86666667e-03j
8.13333333e-03+3.26666667e-03j 1.35000000e-03+3.33333333e-05j
4.68333333e-03+2.36666667e-03j 2.53333333e-03+2.66666667e-04j
3.63333333e-03+8.66666667e-04j 4.93333333e-03+1.86666667e-03j
6.43333333e-03+3.26666667e-03j 8.13333333e-03+5.06666667e-03j
8.83333333e-04+3.33333333e-05j 4.68333333e-03+3.43333333e-03j
1.63333333e-03+2.66666667e-04j 2.53333333e-03+8.66666667e-04j
3.63333333e-03+1.86666667e-03j 4.93333333e-03+3.26666667e-03j
6.43333333e-03+5.06666667e-03j 8.13333333e-03+7.26666667e-03j
5.16666667e-04+3.33333333e-05j 4.68333333e-03+4.70000000e-03j
9.33333333e-04+2.66666667e-04j 1.63333333e-03+8.66666667e-04j
2.53333333e-03+1.86666667e-03j 3.63333333e-03+3.26666667e-03j
4.93333333e-03+5.06666667e-03j 6.43333333e-03+7.26666667e-03j
8.13333333e-03+9.86666667e-03j 2.50000000e-04+3.33333333e-05j
4.68333333e-03+6.16666667e-03j 4.33333333e-04+2.66666667e-04j
9.33333333e-04+8.66666667e-04j 1.63333333e-03+1.86666667e-03j
2.53333333e-03+3.26666667e-03j 3.63333333e-03+5.06666667e-03j
4.93333333e-03+7.26666667e-03j 6.43333333e-03+9.86666667e-03j
8.13333333e-03+1.28666667e-02j 8.33333333e-05+3.33333333e-05j
4.68333333e-03+7.83333333e-03j 1.33333333e-04+2.66666667e-04j
4.33333333e-04+8.66666667e-04j 9.33333333e-04+1.86666667e-03j
1.63333333e-03+3.26666667e-03j 2.53333333e-03+5.06666667e-03j
3.63333333e-03+7.26666667e-03j 4.93333333e-03+9.86666667e-03j
6.43333333e-03+1.28666667e-02j 8.13333333e-03+1.62666667e-02j
1.25000000e-05+2.50000000e-05j 3.09583333e-03+6.19166667e-03j
1.66666667e-05+1.66666667e-04j 1.33333333e-04+8.66666667e-04j
4.33333333e-04+1.86666667e-03j 9.33333333e-04+3.26666667e-03j
1.63333333e-03+5.06666667e-03j 2.53333333e-03+7.26666667e-03j
3.63333333e-03+9.86666667e-03j 4.93333333e-03+1.28666667e-02j
6.43333333e-03+1.62666667e-02j 3.91666667e-03+9.36666667e-03j
1.66666667e-05+5.00000000e-04j 1.33333333e-04+1.86666667e-03j
4.33333333e-04+3.26666667e-03j 9.33333333e-04+5.06666667e-03j
1.63333333e-03+7.26666667e-03j 2.53333333e-03+9.86666667e-03j
3.63333333e-03+1.28666667e-02j 4.93333333e-03+1.62666667e-02j
3.08333333e-03+9.36666667e-03j 1.66666667e-05+1.03333333e-03j
1.33333333e-04+3.26666667e-03j 4.33333333e-04+5.06666667e-03j
9.33333333e-04+7.26666667e-03j 1.63333333e-03+9.86666667e-03j
2.53333333e-03+1.28666667e-02j 3.63333333e-03+1.62666667e-02j
2.35000000e-03+9.36666667e-03j 1.66666667e-05+1.76666667e-03j
1.33333333e-04+5.06666667e-03j 4.33333333e-04+7.26666667e-03j
9.33333333e-04+9.86666667e-03j 1.63333333e-03+1.28666667e-02j
2.53333333e-03+1.62666667e-02j 1.71666667e-03+9.36666667e-03j
1.66666667e-05+2.70000000e-03j 1.33333333e-04+7.26666667e-03j
4.33333333e-04+9.86666667e-03j 9.33333333e-04+1.28666667e-02j
1.63333333e-03+1.62666667e-02j 1.18333333e-03+9.36666667e-03j
1.66666667e-05+3.83333333e-03j 1.33333333e-04+9.86666667e-03j
4.33333333e-04+1.28666667e-02j 9.33333333e-04+1.62666667e-02j
7.50000000e-04+9.36666667e-03j 1.66666667e-05+5.16666667e-03j
1.33333333e-04+1.28666667e-02j 4.33333333e-04+1.62666667e-02j
4.16666667e-04+9.36666667e-03j 1.66666667e-05+6.70000000e-03j
1.33333333e-04+1.62666667e-02j 1.83333333e-04+9.36666667e-03j
1.66666667e-05+8.43333333e-03j 5.00000000e-05+9.36666667e-03j
4.16666667e-06+3.17500000e-03j]
定义Dirichlet边界条件，建立并求解变分问题。
求解变分问题#
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
bc = dolfinx.fem.dirichletbc(u_c, boundary_dofs)
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
uh = problem.solve()
计算 \(L^2\) 误差和最大误差
误差计算#
x = ufl.SpatialCoordinate(mesh)
u_ex = 0.5 * x[0]**2 + 1j*x[1]**2
L2_error = dolfinx.fem.form(ufl.dot(uh-u_ex, uh-u_ex) * ufl.dx(metadata={"quadrature_degree": 5}))
local_error = dolfinx.fem.assemble_scalar(L2_error)
global_error = np.sqrt(mesh.comm.allreduce(local_error, op=MPI.SUM))
max_error = mesh.comm.allreduce(np.max(np.abs(u_c.x.array-uh.x.array)))
print(global_error, max_error)
(0.0007865435216227265+0.0017660156338113752j) 3.553078358632328e-06
绘图#
分别画出解的实部和虚部。
import pyvista
pyvista.start_xvfb()
p_mesh = pyvista.UnstructuredGrid(*dolfinx.plot.create_vtk_mesh(mesh, mesh.topology.dim))
pyvista_cells, cell_types, geometry = dolfinx.plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u_real"] = uh.x.array.real
grid.point_data["u_imag"] = uh.x.array.imag
_ = grid.set_active_scalars("u_real")
p_real = pyvista.Plotter()
p_real.add_text("uh real", position="upper_edge", font_size=14, color="black")
p_real.add_mesh(grid, show_edges=True)
p_real.view_xy()
if not pyvista.OFF_SCREEN:
p_real.show()
grid.set_active_scalars("u_imag")
p_imag = pyvista.Plotter()
p_imag.add_text("uh imag", position="upper_edge", font_size=14, color="black")
p_imag.add_mesh(grid, show_edges=True)
p_imag.view_xy()
if not pyvista.OFF_SCREEN:
p_imag.show()