You can use Feynman-Kac formula:
$$ u(x) = -\mathbb{E}_x\int_0^T f(X_t)dt$$where $X_t$ is a Brownian process starting at $x$ and $T$ the stopping time $T=\inf\{t\geq 0:X(t)\in \partial D\}$. So $u(x)$ is essentially given by the expected time for the process to get out of $D$. For example if $\epsilon \leq f\leq M$ we get $$ \epsilon\mathbb{E}_x(T)\leq -u(x) \leq M \mathbb{E}_x(T).$$
Consider the case $\partial D$ smooth around a point $x_0$. As a simplification we suppose that locally $D=\mathbb{R}_+\times\mathbb{R}^{n-1}$, $x_0 = (0,\cdots,0)$ and $x=(\delta,0,\cdots,0)$ with $\delta>0$ and $X_s = (x+B_s^{(1)},B_s^{(2)},\cdots,B_s^{(n)})$, where $B^{(i)}$ are iid Brownian motion. In that case it is easy to gives the escape time of $X$: we have $$\mathbb{P}_x(T>t)=\mathbb{P}_x(\inf_{0\leq s\leq t} B_s^{(1)} > - \delta) = 1-2\mathbb{P}_x(B_t^{(1)}<-\delta) \approx \frac{2\delta}{\sqrt{2\pi t}} $$ (see reflection) and then
$$ \mathbb{E}_x(T)= \mathbb{E}_x(\int_0^\infty 1_{t< T}dt)=\mathbb{E}_x(\int_0^1 1_{t< T}dt)+\mathbb{E}_x(1_{T>1}\int_1^\infty 1_{t\leq T}dt) \\ = \int_0^1 \mathbb{P}_x(t< T)dt+\mathbb{P}_x(T>1) \mathbb{E}_x\left( \int_1^\infty 1_{t\leq T}dt|T>1\right)
\approx \delta C$$
for some $C>0$.
If the boundary is not regular at $x_0$, it can be possible for the process to avoid $\partial D$ such that $\mathbb{P}(T>1)$ is bounded away from $0$ (or decay very slowly) as $x\rightarrow x_0$ and then $u(x)$ is no more comparable with $d(x,\partial D)$. This lead to a nice physical phenomena call " Effet de pointe" that occures in Lightning rod.