如何修复FiPy中矩形网格下多元高斯初始条件失效问题

我尝试在FiPy数值模拟中使用多元高斯初始条件，当前代码在方形网格下可生成正常的高斯分布效果，但应用于矩形网格时效果异常。代码如下：

from fipy import CellVariable, Grid2D, Viewer
from scipy.stats import multivariate_normal
import numpy as np
import matplotlib.pyplot as plt
plt.close('all')
# Define the grid and cell variable
nx = 40
ny = 100
dx = 1.0
dy = 1.90
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = CellVariable(name="phi", mesh=mesh)

# Set the Gaussian initial condition
mean = [nx * dx / 2, ny * dy / 2]  # Center of the Gaussian surface
covariance = [[10, 0], [0, 5]]  # Covariance matrix

# Generate coordinates for the grid
X, Y = mesh.cellCenters[0], mesh.cellCenters[1]

# Evaluate the Gaussian surface
gaussian_surface = multivariate_normal(mean=mean, cov=covariance)
Z = np.zeros_like(X)
Xm=X.value.reshape([nx,ny])
Ym=Y.value.reshape([nx,ny])

Z = gaussian_surface.pdf(np.column_stack((X.value.flat, Y.value.flat)))

# Assign the Gaussian surface to the cell variable
phi.setValue(Z)

plt.pcolor(Xm,Ym,phi.value.reshape((nx, ny)), cmap='plasma')

plt.colorbar(label='phi')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Gaussian Initial Condition')
plt.show()

修复方案 #

矩形网格下效果异常的核心原因有两个，对应两种修复方向：

1. 可视化坐标处理不当 #

plt.pcolor 默认需要传入网格节点坐标，但原代码使用的是单元中心坐标，且未匹配物理尺度的坐标轴比例，导致显示变形。

修复方法：

• 使用FiPy自带Viewer（推荐，自动处理网格坐标）：
  替换原plt相关代码为：

  viewer = Viewer(vars=phi, cmap='plasma')
  viewer.plot()
  plt.title('Gaussian Initial Condition')
  plt.show()
  

• 调整matplotlib可视化逻辑：
  传入节点坐标并设置等物理尺度坐标轴：

  # 获取网格节点坐标（节点数=单元数+1）
  X_nodes = mesh.faceCenters[0].value.reshape(nx+1, ny+1)
  Y_nodes = mesh.faceCenters[1].value.reshape(nx+1, ny+1)
  # 绘制并启用等物理尺度显示
  plt.pcolor(X_nodes, Y_nodes, phi.value.reshape(nx, ny), cmap='plasma')
  plt.colorbar(label='phi')
  plt.xlabel('X')
  plt.ylabel('Y')
  plt.title('Gaussian Initial Condition')
  plt.axis('equal')
  plt.show()
  

2. 协方差未匹配网格尺度 #

若期望高斯分布的形状相对于网格单元数保持一致（比如方形网格下x方向10个单元、y方向5个单元的分布宽度），原代码的协方差是基于物理坐标的，矩形网格下dx≠dy会导致物理空间的分布比例失衡。

修复方法：

将协方差矩阵的元素乘以对应方向网格间距的平方，把单元数尺度转换为物理坐标尺度：

# 修改协方差矩阵定义
covariance = [[10 * dx**2, 0], [0, 5 * dy**2]]

完整修正代码（保持单元数比例+FiPy Viewer可视化） #

from fipy import CellVariable, Grid2D, Viewer
from scipy.stats import multivariate_normal
import numpy as np
import matplotlib.pyplot as plt
plt.close('all')

# Define the grid and cell variable
nx = 40
ny = 100
dx = 1.0
dy = 1.90
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = CellVariable(name="phi", mesh=mesh)

# Set the Gaussian initial condition
mean = [nx * dx / 2, ny * dy / 2]
# 协方差转换为物理尺度（匹配单元数比例）
covariance = [[10 * dx**2, 0], [0, 5 * dy**2]]

# Generate coordinates for the grid
X, Y = mesh.cellCenters[0], mesh.cellCenters[1]

# Evaluate the Gaussian surface
gaussian_surface = multivariate_normal(mean=mean, cov=covariance)
Z = gaussian_surface.pdf(np.column_stack((X.value.flat, Y.value.flat)))

# Assign the Gaussian surface to the cell variable
phi.setValue(Z)

# 使用FiPy Viewer可视化
viewer = Viewer(vars=phi, cmap='plasma')
viewer.plot()
plt.title('Gaussian Initial Condition')
plt.show()

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