defining a multivariate gaussian initial condition in FiPy

I'm trying to use a multivariate gaussian initial condition for a Fipy integration.

I'm currently using the following code:

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()

The code I have works well for square grids:

But It does not work well for rectangular ones:

How can I fix it?

The OP's code gives the following warning:

/tmp/ipykernel_9302/2008401024.py:32: UserWarning: The input coordinates to 
pcolor are interpreted as cell centers, but are not monotonically 
increasing or decreasing. This may lead to incorrectly calculated 
cell edges, in which case, please supply explicit cell edges to 
pcolor.  plt.pcolor(Xm,Ym,phi.value.reshape((nx, ny)), cmap='plasma')

This hints at a reshaping problem, due to row-major matrix storage being treated as column-major (or vice-versa). So, to fix it, we can simply reverse the (nx, ny) in the final steps using nx, ny = ny, nx. As expected, the OP's original code worked fine when nx=ny. The corrected version is shown below for completeness:

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)

nx, ny = ny, nx

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()

Output: