使用 fipy 在坐标上定义的网格进行生物扩散的表示问题

Question

我想使用 fipy 通过生物扩散来模拟狼群动力学。我使用了来自荷兰的数据，因此我的网格是根据坐标定义的，不再从零开始。初始条件给了我预期的结果（其他地方 x0、y0 和 0 处的最大密度），但后续快照显示整个网格的密度为零。

我创建了一个初始 2D 模型，通过考虑从 0 延伸到 Lx 的域，该模型可以完美地工作。该模型根据初始条件向我展示了某个位置的狼群，然后该狼群根据其他狼群（本身是固定的）引起的排斥力移动。因此，我获得了几张快照，显示狼的密度正在移动。现在我已经根据坐标更改了网格和轴，初始条件是正确的，但后续快照显示整个区域的密度为零，而我期望仅通过更改坐标获得与基本模型类似的结果.

from fipy import CellVariable, Grid2D, ExponentialConvectionTerm, TransientTerm, DiffusionTerm
from fipy.tools import numerix
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

PB = plt.imread("Pays-Bas.png")

# Model parameters
xs, xe = 684000, 705878      # X coordinates of the rectangle surrounding the area of interest [m] (WGS84 - UTM31) 
ys, ye = 5768000, 5779500    # Y coordinates [m]  (WGS84 - UTM31)
Lx = xe - xs
Ly = ye - ys
s0 = 100
x0, y0 = 693822.958, 5775254.324   
Smax = 1                          

# Other packs
xPack1, yPack1 = 688288.499, 5779104.309
xPack2, yPack2 = 684451.539, 5768900.565
xPack3, yPack3 = 705641.974, 5771640.434
repulsion_strength = 50
repulsion_range = np.sqrt(180/np.pi)*1000

# Numerical parameters
Nx, Ny = 550, 300
dx, dy = Lx/Nx, Ly/Ny
T = 7*24*3600.  # Simulation sur 7 jours
dt = T/200
Np = 7
Nt = Np*np.ceil(T/(Np*dt)).astype('int')

# Define the grid/mesh
mesh = Grid2D(nx=Nx, ny=Ny, dx=dx, dy=dy)
x, y = mesh.cellCenters[0], mesh.cellCenters[1]

xg, yg = np.meshgrid(np.linspace(xs, xe, Nx+1), np.linspace(ys, ye, Ny+1))
xd, yd = np.meshgrid(np.linspace(xs + dx/2, xe - dx/2, Nx),
                     np.linspace(ys + dy/2, ye - dy/2, Ny))

# Define the model variable and set the initial and boundary conditions
def ic(x, y):
    return Smax*numerix.exp(-((x - x0)**2 + (y - y0)**2)/(2*s0**2))

S = CellVariable(name="numerical solution", mesh=mesh, hasOld=True, value=ic(x, y))
S.faceGrad.constrain(0, mesh.facesLeft)
S.faceGrad.constrain(0, mesh.facesRight)
S.faceGrad.constrain(0, mesh.facesTop)
S.faceGrad.constrain(0, mesh.facesBottom)

# Visualisation de la condition initiale
fig_init, ax_init = plt.subplots(1, 1)
ax_init.imshow(PB, extent=[xs, xe, ys, ye])
c_init = ax_init.pcolormesh(xg, yg, ic(xd,yd), cmap='jet')
plt.colorbar(c_init)
plt.title("Condition Initiale")
plt.scatter([xPack1, xPack2, xPack3], [yPack1, yPack2, yPack3], marker='o', color='red', label='Meute adverse')
plt.scatter(xPhi, yPhi, marker='o', color='green', label='Target')
plt.imshow(PB, zorder=1, alpha=0.5, extent=[xs,xe,ys,ye])
plt.show()

# Define and then solve the equation
def phi(x, y):
    repulsionPack1 = repulsion_strength*numerix.exp(-((x - xPack1)**2 + (y - yPack1)**2)/(2*repulsion_range**2))
    repulsionPack2 = repulsion_strength*numerix.exp(-((x - xPack2)**2 + (y - yPack2)**2)/(2*repulsion_range**2))
    repulsionPack3 = repulsion_strength*numerix.exp(-((x - xPack3)**2 + (y - yPack3)**2)/(2*repulsion_range**2))
    return repulsionPack1 + repulsionPack2 + repulsionPack3

phi_var = CellVariable(name="phi", mesh=mesh, value=phi(x, y))
eq = TransientTerm() == DiffusionTerm(coeff=alpha + beta*S) \
    + ExponentialConvectionTerm(coeff=phi_var.faceGrad)

# Matrix with Np solution snapshots
my_sol = np.zeros((Np, Nx*Ny))
my_sol[0, :] = S

k = 1

for step in np.arange(1, Nt):
    S.updateOld()

    res = 1.
    while res > 1.e-6:
        res = eq.sweep(var=S, dt=dt)

    if np.mod(step, Nt/Np) == 0:
        print(f"Snapshot {k} à t = {step*dt/(24*3600)} jours")
        my_sol[k, :] = S
        k += 1

print(my_sol)

# Plot 2D
fig, ax = plt.subplots(1,1)
ax.imshow(PB, extent=[xs, xe, ys, ye])

for i in np.arange(Np):
    sol = my_sol[i, :].reshape((Nx, Ny))
    sol = griddata((xd.ravel(), yd.ravel()), my_sol[i, :], (xg, yg), method='nearest')
    c = ax.pcolormesh(xg, yg, sol, cmap='jet')
    plt.colorbar(c)
    plt.title(f"Snapshot {i + 1}")
    plt.legend(fontsize='small',ncol=1)
    plt.scatter([xPack1, xPack2, xPack3],[yPack1, yPack2, yPack3], marker='o', color='red', label='Meute adverse')
    plt.imshow(PB, zorder=1, alpha=0.5, extent=[xs,xe,ys,ye])
    plt.show()

Answer 1

我认为你只需要偏移网格以与你的数据一致：

mesh = Grid2D(nx=Nx, ny=Ny, dx=dx, dy=dy) + [[xs], [ys]]