Another approach that confirms the above extremely close approximation is to introduce
$$S = \log P = \log \prod_{n\ge 1} \left(1-\frac{1}{2^n}\right)
= \sum_{n\ge 1} \log \left(1-\frac{1}{2^n}\right)$$
and observe that this sum is harmonic and may be evaluated by inverting its Mellin transform. To do this introduce
$$S(x) = \sum_{n\ge 1} \log \left(1-\frac{1}{2^{nx}}\right)$$
Recall the harmonic sum identity
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad
g(x) = \log\left(1-\frac{1}{2^x}\right).$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is
$$\int_0^\infty \log\left(1-\frac{1}{2^x}\right) x^{s-1} dx.$$
The function $g(x)$ is well-behaved near zero where it is on the order of $\log x$ and vanishes faster than any polynomial at infinity.
To calculate the Mellin transform start with
$$\int_0^\infty \log\left(1-\frac{1}{2^x}\right) x^{s-1} dx
= - \int_0^\infty \sum_{q\ge 1} \frac{2^{-qx}}{q} x^{s-1} dx
= - \sum_{q\ge 1} \frac{1}{q} \int_0^\infty 2^{-qx} x^{s-1} dx.$$
Observe that
$$\int_0^\infty 2^{-qx} x^{s-1} dx =
\frac{1}{(q \log 2)^s} \Gamma(s)$$
by a straightforward substitution that turns the integral into a gamma function integral.
This yields
$$g^*(s) = - \sum_{q\ge 1} \frac{1}{q} \frac{1}{(q \log 2)^s} \Gamma(s)
= -\frac{1}{(\log 2)^s} \Gamma(s) \sum_{q\ge 1} \frac{1}{q^{s+1}}
= -\frac{1}{(\log 2)^s} \Gamma(s) \zeta(s+1).$$
By the harmonic sum identity we now have that the Mellin transform $Q(s)$ of $S(x)$ is given by $$Q(s) = -\frac{1}{(\log 2)^s} \Gamma(s) \zeta(s) \zeta(s+1)$$
with the Mellin inversion integral being
$$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$
which we evaluate by shifting it to the left for an expansion about zero.
Fortunately the two zeta function terms with their trivial zeros combine to cancel the poles of the gamma function and we are left with just three poles and residues.
We have
$$\mathrm{Res}(Q(s)/x^s; s=1) = -\frac{\pi^2}{6x\log 2}.$$
Furthermore
$$\mathrm{Res}(Q(s)/x^s; s=0) =
\frac{1}{2}\left(\log\frac{2\pi}{\log 2}-\log x\right)$$
and finally
$$\mathrm{Res}(Q(s)/x^s; s=-1) = \frac{1}{24} x \log 2.$$
Putting $x=1$ we obtain the following approximation for $S(1):$
$$S(1)\approx \frac{1}{24} \log 2 + \frac{1}{2} \log\frac{2\pi}{\log 2}
- \frac{\pi^2}{6\log 2}.$$
This is $$-1.2420620948124149457978452979784311762117047031228$$
while the exact value is
$$-1.2420620948124149457978454818946296689734039782504$$
so this approximation is good to an amazing $25$ digits.
This gives for $P$ the approximation
$$P \approx 2^{1/24} \sqrt{\frac{2\pi}{\log 2}}
\exp\left(- \frac{\pi^2}{6\log 2}\right).$$
This is also good to $25$ digits, confirming the observation from the other poster above.