求解无电荷2D泊松方程时，如何在FiPy中计算极化矢量散度？

要在FiPy中求解无电荷二维泊松方程 $\nabla·\epsilon_0 \nabla\phi = \nabla·\vec{P}$，关键是正确计算矢量场$\vec{P}$的散度，具体操作步骤如下：

• 先确保$\vec{P}$被定义为rank=1的CellVariable（对应二维矢量场）：

  from fipy import Grid2D, CellVariable
  
  mesh = Grid2D(nx=50, ny=50)  # 根据需求自定义网格规模
  P = CellVariable(mesh=mesh, rank=1)  # rank=1标记该变量为矢量场
  # 可在此处为P赋值，例如 P.setValue([P_x, P_y])
  

• 计算$\vec{P}$的散度需先将单元矢量插值为面变量（FiPy的散度计算基于面通量离散实现），再调用divergence()函数：

  from fipy import divergence
  
  # 将单元矢量P插值到网格面
  P_face = P.faceValue
  # 计算矢量场的散度
  div_P = divergence(P_face)
  

• 最后构建并求解电压方程：

  from fipy import DiffusionTerm
  
  epsilon0 = 8.854e-12
  phi = CellVariable(mesh=mesh)
  # 构建对应泊松方程的FiPy形式
  eq = DiffusionTerm(coeff=epsilon0, var=phi) == div_P
  # 根据问题需求添加边界条件（如Dirichlet或Neumann条件）
  # 执行方程求解
  eq.solve()
  

关键说明 #

• 你之前尝试的u.divergence无效，是因为FiPy中divergence()是独立函数而非矢量变量的内置方法，且它仅作用于面矢量变量，无法直接处理单元矢量。
• 手册中提到的var.faceGrad.divergence本质是对标量场的梯度（面矢量）求散度，和矢量场直接求散度的逻辑一致：先得到面矢量，再计算散度。

内容的提问来源于stack exchange，提问作者Kritika Khanal