Thanks for the reply. I still could not solve the problem. The model looks good and the code is all right. Perhaps with a long time computation, the outcome is okay, but it is not efficient.
I encountered this problem a few weeks ago when I was reading the reference “Flexible modelling of spatial variation in agricultural field trials with the R package INLA” and tried to use INLA to fit the model with AR1\otimes AR1 spatial residual covariance. I used brms in the first place and found that it is not capable to do so. Then I coded up in rstan myself and am facing this problem.