Please search the forum for similar questions, for instance Fenicsx: use an array of data for boundary conditions (with a couple of links) or Mapping 2D numpy array into dolfinx function
I wouldn’t trust
# applying volume source
ids = msh.geometry.input_global_indices
dim = msh.topology.dim
data = np.loadtxt("source.raw")
source.x.array[:] = data[ids, 3:3+dim].flatten()
as there is nothing that guarantees you that the order of the points in the file is the same at the locations of the DOFs.