def build_derivatives_2d(NAS, NVS, ds, dv): N = NAS * NVS DS = lil_matrix((N, N)) D2S = lil_matrix((N, N)) DV = lil_matrix((N, N)) D2V = lil_matrix((N, N)) DSV = lil_matrix((N, N)) def idx(i, j): # Flattening return i + j * NAS for j in range(NVS): for i in range(NAS): k = idx(i, j) # First derivative S if 1 <= i < NAS - 1: DS[k, idx(i - 1, j)] = -0.5 / ds DS[k, idx(i + 1, j)] = 0.5 / ds elif i == 0: DS[k, idx(i, j)] = -1.0 / ds DS[k, idx(i + 1, j)] = 1.0 / ds elif i == NAS - 1: DS[k, idx(i - 1, j)] = -1.0 / ds DS[k, idx(i, j)] = 1.0 / ds # Second derivative S if 1 <= i < NAS - 1: D2S[k, idx(i - 1, j)] = 1.0 / ds**2 D2S[k, idx(i, j)] = -2.0 / ds**2 D2S[k, idx(i + 1, j)] = 1.0 / ds**2 elif i == 0: D2S[k, idx(i, j)] = 1.0 / ds**2 D2S[k, idx(i + 1, j)] = -2.0 / ds**2 D2S[k, idx(i + 2, j)] = 1.0 / ds**2 elif i == NAS - 1: D2S[k, idx(i - 2, j)] = 1.0 / ds**2 D2S[k, idx(i - 1, j)] = -2.0 / ds**2 D2S[k, idx(i, j)] = 1.0 / ds**2 # First derivative v if 1 <= j < NVS - 1: DV[k, idx(i, j - 1)] = -0.5 / dv DV[k, idx(i, j + 1)] = 0.5 / dv elif j == 0: DV[k, idx(i, j)] = -1.0 / dv DV[k, idx(i, j + 1)] = 1.0 / dv elif j == NVS - 1: DV[k, idx(i, j - 1)] = -1.0 / dv DV[k, idx(i, j)] = 1.0 / dv # Second derivative if 1 <= j < NVS - 1: D2V[k, idx(i, j - 1)] = 1.0 / dv**2 D2V[k, idx(i, j)] = -2.0 / dv**2 D2V[k, idx(i, j + 1)] = 1.0 / dv**2 elif j == 0: D2V[k, idx(i, j)] = 1.0 / dv**2 D2V[k, idx(i, j + 1)] = -2.0 / dv**2 D2V[k, idx(i, j + 2)] = 1.0 / dv**2 elif j == NVS - 1: D2V[k, idx(i, j - 2)] = 1.0 / dv**2 D2V[k, idx(i, j - 1)] = -2.0 / dv**2 D2V[k, idx(i, j)] = 1.0 / dv**2 # Cross derivative (central differences) if (1 <= i < NAS - 1) and (1 <= j < NVS - 1): DSV[k, idx(i - 1, j - 1)] = 1.0 / (4 * ds * dv) DSV[k, idx(i + 1, j + 1)] = 1.0 / (4 * ds * dv) DSV[k, idx(i - 1, j + 1)] = -1.0 / (4 * ds * dv) DSV[k, idx(i + 1, j - 1)] = -1.0 / (4 * ds * dv) return DS.tocsr(), D2S.tocsr(), DV.tocsr(), D2V.tocsr(), DSV.tocsr()

Apologies, but something went wrong on our end.