import numpy as np

import matplotlib.pyplot as plt

# Grid and physical parameters

Nx, Ny = 50, 50

Lx, Ly = 1.0, 1.0

dx, dy = Lx / (Nx - 1), Ly / (Ny - 1)

Ra = 1e5

Pr = 0.71

nu = 1.0 / np.sqrt(Ra * Pr)

alpha = nu / Pr

dt = 0.001

# Meshgrid

x = np.linspace(0, Lx, Nx)

y = np.linspace(0, Ly, Ny)

X, Y = np.meshgrid(x, y, indexing='ij')

# Triangular mask: only compute where x + y

Andrew Ssenabulya, Thank you for your response.