In general this will not have a solution. Note that your equation is equivalent to the following:
\[ \Delta_g (u) = \sum_{i} \frac{1}{a} \frac{\partial}{\partial \theta_i}( a \frac{\partial u}{\partial \theta_i} ) = e^{i\theta_1} \]
where by definition $a = e^b$. This is the Laplace equation (see the Wikipedia article "Laplace-Beltrami operator") with respect to the metric
\[ g_{ij} = e^{2b} \delta_{ij} \]
Hence it only has a solution if the average value of $e^{i\theta_1}$ with respect to the volume form of $g$ is zero. The reason for this is the integration by parts formula
\[ \int_{[0,2\pi]^d} (\Delta_g u)\cdot v dvol_g = \int_{[0,2\pi]^d} u \cdot (\Delta_g v) dvol_g \]
which when you plug in $v = 1$ shows that the average value of $\Delta_g u $ must be zero. The formula can be proved in much the same way as for the standard Laplacian. Thus a necessary (and sufficient) condition for a solution to exist is:
\[ \int_{[0,2\pi]^d} e^{b+ i \theta_1} d\theta_1 \cdots d\theta_n = 0\]
In general, even if a solution exists I wouldn't expect a particularly nice description of the Fourier coefficients. You might be able to get estimates of them, particularly if the function $b$ is very small, in which case the equation will closely resemble the standard Laplacian on the torus. If $b$ has very few terms you might try taking the Fourier transform of your equation, in which case the second term will become a convolution and you may be able to solve inductively for the Fourier coefficients if you're lucky.
