It may be of interest to note that the sum
$$ g(x) = \sum_{k=1}^\infty \log\tanh (kx)$$
is harmonic and may be evaluated using Mellin transforms, yielding an asymptotic expansion about zero.
The Mellin transform $f^*(s)$ of the base function
$$ f(x) = \log\tanh x$$
may be computed as follows
\begin{align}
f^*(s) & = \mathfrak{M}(f(x); s) =
\int_0^\infty \log \frac{e^x-e^{-x}}{e^x+e^{-x}} x^{s-1} \, dx
= \int_0^\infty \log \left(1 - 2\frac{e^{-x}}{e^x+e^{-x}}\right) x^{s-1} dx \\[6pt]
& = \int_0^\infty \log \left(1 - 2\frac{e^{-2x}}{1+e^{-2x}}\right) x^{s-1} dx
= - \int_0^\infty \sum_{q\ge 1} \frac{1}{q} 2^q e^{-2qx}\left(\frac{1}{1+e^{-2x}}\right)^q x^{s-1} dx \\[6pt]
& = -\sum_{q\ge 1} \frac{1}{q} 2^q \int_0^\infty e^{-2qx}
\sum_{m\ge 0} (-1)^m \binom{m+q-1}{m} e^{-2mx} x^{s-1} dx \\[6pt]
& = -\sum_{q\ge 1} \frac{1}{q} 2^q \sum_{m\ge 0} (-1)^m \binom{m+q-1}{m}
\int_0^\infty e^{-2qx} e^{-2mx} x^{s-1} dx \\[6pt]
& =-\Gamma(s)\sum_{q\ge 1} \frac{1}{q} 2^q \sum_{m\ge 0} (-1)^m \binom{m+q-1}{m}
\frac{1}{(2m+2q)^s} \\[6pt]
& = -\frac{\Gamma(s)}{2^s} \sum_{q\ge 1} \frac{1}{q} 2^q \sum_{m\ge 0} (-1)^m \binom{m+q-1}{m} \frac{1}{(m+q)^s}.
\end{align}
To complete this calculation, ask about the coefficient of
$$\frac{1}{n^s} = \frac{1}{(m+q)^s}.$$
It is given by
$$\sum_{m=0}^{n-1} \frac{1}{n-m} 2^{n-m} (-1)^m \binom{n-1}{m} =
- \frac{1}{n} (-1)^n +
\frac{1}{n} \sum_{m=0}^n 2^{n-m} (-1)^m \binom{n}{m} =
\frac{1}{n} \left(1-(-1)^n\right).$$
It follows that the double sum is
$$\sum_{n\ge 1} \frac{1-(-1)^n}{n^{s+1}}=
2 \sum_{m\ge 0} \frac{1}{(2m+1)^{s+1}} =
2 \zeta(s+1) \left(1 - \frac{1}{2^{s+1}}\right) =
\zeta(s+1) \left(2 - \frac{1}{2^s}\right).$$
This gives the following for $f^*(s):$
$$ f^*(s) = -\frac{\Gamma(s)}{2^s} \zeta(s+1) \left(2 - \frac{1}{2^s}\right).$$
Now introduce $g(x)$, the harmonic sum we are trying to calculate, so that
$$ g(x) = \sum_{k=1}^\infty f(kx).$$
The Mellin transform of $g(x)$ is then given by
$$ g^*(s) = \mathfrak{M}(g(x); s) =
-\frac{\Gamma(s)}{2^s} \zeta(s) \zeta(s+1) \left(2 - \frac{1}{2^s}\right).$$
The zeros of the two zeta function terms cancel the poles of the gamma function, so that inverting $g^*(s)$ we only have two terms that contribute, namely
$$\operatorname{Res}(g^*(s) x^{-s}; s=1) = -1/8\,{\frac {{\pi }^{2}}{x}}$$ and
$$\operatorname{Res}(g^*(s) x^{-s}; s=0) =
1/2\,\log \left( 2\,\pi \right) -1/2\,\log \left( x \right).$$
This yields that in a neighborhood of zero
$$ g(x) \sim 1/2\,\log \left( 2\,\pi \right) -1/2\,\log \left( x \right) -1/8\,{\frac {{\pi }^{2}}{x}}$$
Setting $x=\frac{\pi}{2}$, we obtain that
$$g\left(\frac{\pi}{2}\right) \sim \log 2 - \frac{\pi}{4},$$
which produces only three good digits. On the other hand, for e.g. $x=\frac{1}{4}$, we get
$$g\left(\frac{1}{4}\right) \sim 3/2\,\log \left( 2 \right) +1/2\,\log \left( \pi \right) -1/2\,{\pi }^{2} \sim -3.322716487,$$
which has nine good digits.
For e.g. $x=\frac{\pi}{16}$, we get
$$g\left(\frac{\pi}{16}\right) \sim 5/2\,\log \left( 2 \right) -2\,\pi
\sim -4.5503173557797232034$$
which has 20 good digits.
It seems quite intriguing to ask whether this expansion can also be derived directly from properties of the Jacobi theta function without using Mellin transforms.