Solution of Burgers' equation in preshock region According to Hamilton's & Blackstock's Nonlinear acoustics (Section 4.5.4) the solution of Burgers' equation of the form:
$$
\frac{\partial P}{\partial \sigma} - \frac{1}{\Gamma}\frac{\partial^2 P}{\partial \theta^2}=P\frac{\partial P}{\partial \theta}
$$
for source-problem of monofrequency radiation at $x=0$ (amplitude $p_0$, angular frequency $\omega$)
$$
p = p_0 \frac{4\Gamma^{-1}\sum_{n=1}^\infty (-1)^{n+1}nI_n(\frac{1}{2}\Gamma)e^{-n^2\alpha x}\sin n\omega\tau}{I_0(\frac{1}{2}\Gamma)+2\sum_{n=1}^\infty (-1)^{n}I_n(\frac{1}{2}\Gamma)e^{-n^2\alpha x}\cos n\omega\tau}
$$
where


*

*$P = \frac{p}{p_0}$

*$\sigma = \frac{x}{\bar{x}}$, $\ \ \bar{x}$ being shock formation distance

*$\theta = \omega \tau$,  $ \ \ \tau$ being retarded time

*$\Gamma = \frac{\ell_a}{\bar{x}}$, $\ \ \ell_a$ being absorption length, $\ \ \ell_a = \frac{1}{\alpha}$, where $\alpha$ is the absorption coefficient

*$I_n(x)$ modified Bessel functions


This solution should work for all distances (i.e. for all stages of shock growth, shock decay etc.). It works fine for the dissipation of shock in the sawtooth region ($\sigma > 3$) but I am not successful in the preshock region.
The expression is independent of $\sigma$ and therefore there is no higher harmonics growth etc., only decay of them. Does that mean that this solution is applicable only in the state when all the harmonics are present (already grown) and it maps only their decay according to mode number $n$?
 A: 
The expression is independent of $\sigma$ and therefore there is no higher harmonics growth etc., only decay of them.

I think your exponents have $\sigma$ in them through the $x$ terms, e.g.,
$$
e^{-n^{2} \alpha \ x} = e^{-n^{2} \alpha \ \bar{x} \ \sigma} = e^{-n^{2} \frac{\sigma}{\Gamma}}
$$
One could also find $\sigma$ from $\Gamma = \tfrac{\sigma \ l_{a}}{x}$.
As a side note, a shock wave is a discontinuity, which means to construct it one must use an infinite number of Fourier terms.  In other words, the sawtooth description you describe already has an infinite number of harmonics for it must to create a discontinuity.

It works fine for the dissipation of shock in the sawtooth region ($\sigma > 3$) but I am not successful in the preshock region.

I am not sure I follow this argument since the $\partial_{\sigma}$ and $\partial_{\theta}$ terms should go to zero in the upstream(preshocked) region, which satisfies the Burgers' equation perfectly well.  Note that $\partial_{\lambda}$ is just short hand for $\tfrac{\partial}{\partial \lambda}$, not a distinction between covariant vs. contravariant derivatives.

Does that mean that this solution is applicable only in the state when all the harmonics are present (already grown) and it maps only their decay according to mode number $n$?

No, I think you missed the $\sigma$ terms as I mentioned above.  However, the Burgers' equation is used to describe a nonlinear wave with a form of diffusion or viscosity as the dissipative term (e.g., $\nu \sim \Gamma^{-1}$ in your equation).  In the region where the shock/discontinuity does not exist, all $\partial_{\lambda}$ terms go to zero so the equation still holds.
Example Result
I used your solution for $P$ in Mathematica and found that it looks perfectly fine to me (see plot output below).
The code I used is given by:
    amp[xb_, aa_, xx_, tt_, nt_] := Module[
      {chi = xb, alpha = aa, x = xx, theta = tt, 
       nn = nt, gam, numer, denom},
      gam = (1/(alpha chi));
      numer = (4/gam) Sum[((-1)^(n+1) n BesselI[n,gam/2] Exp[-n^(2) alpha x] Sin[n theta]),{n,1,nn,1}];
      denom = BesselI[0,gam/2] + 2 Sum[((-1)^(n) BesselI[n,gam/2] Exp[-n^(2) alpha x] Cos[n theta]),{n,1,nn,1}]
      (* output result *)
      N[(numer/denom)]
      ]

I then plotted this using:
    Plot[amp[2, 1, xx, \[Pi]/4, 2000], {xx, -2, 5}, PlotRange -> All, PlotStyle -> {Thick, Red}]

Which gave me the following result:

Side Note
The Burgers' equation is normally written with the temporal and spatial derivatives switched from the form you showed, e.g., given by:
$$
\partial_{\tau} \rho + \rho \partial_{\chi} \rho = \nu \ \partial_{\chi \chi} \rho
$$
where $\rho$ is the amplitude, $\nu$ a constant diffusion or viscosity, $\tau$ a dimensionless temporal variable, and $\chi$ a dimensionless spatial variable.
