Wall function with pressure gradient as described in [42] . Activate with control file command STOKES_BC.

Boundary layer velocity profile:

\begin{eqnarray*} u^+_{P^+=0} &=& \frac{1}{\kappa} \ln \left(1 + \kappa y^+ \right) + \left(C - \frac{1}{\kappa} \ln \kappa \right) \left[1 - \mathrm{e}^{-\frac{y^+}{\eta^+}} - \frac{y^+}{\eta^+} \mathrm{e}^{-by^+}\right] \\ u^+ &=& u^+_{P^+=0} + \left(1 - \mathrm{e}^{-\frac{y^+}{\eta^+}} \right) \cdotp \frac{1}{\kappa} \left[ 2 \left( \sqrt{P^+ y^+ +1} - 1 \right) + \ln 4 - \ln\left(2 + P^+ y^+ + 2 \sqrt{P^+y^++1} \right) \right] \\ \end{eqnarray*}

Matching point of sublayer and log-layer solutions:

\[ \eta^+ + \frac{1}{2} P^+ {\eta^+} ^2 = \frac{1}{\kappa} \ln \eta^+ + C + \frac{1}{\kappa} \left[ 2 \left( \sqrt{P^+ \eta^+ + 1} -1 \right) + \ln 4 - \ln \left( 2 + P^+ \eta^+ + 2 \sqrt{P^+ \eta^+ + 1} \right) \right] \]

For coupled solving of two equations, a two-dimensional Newton Raphson scheme is used. Jacobians are computed per element to simplify matrix and determinant computation.