1= r"""
2 for (int n = 0; n < NumberOfUnknowns; n++) {
3 S[n] = 0.0;
4 }
5
6 if (Q[Shortcuts::h] <= hThreshold) {
7 return;
8 }
9
10 // Compute local slope cF
11 const auto dzx{deltaQ[Shortcuts::z]};
12 const auto dzy{deltaQ[Shortcuts::z + NumberOfUnknowns + NumberOfAuxiliaryVariables]};
13 const auto cF{1.0 / std::sqrt(1.0 + dzx * dzx + dzy * dzy)};
14
15 // Apply friction only if there is some movement
16 const auto momentumSQR{Q[Shortcuts::hu] * Q[Shortcuts::hu] + Q[Shortcuts::hv] * Q[Shortcuts::hv]};
17 const auto momentum{std::sqrt(momentumSQR)};
18 if (tarch::la::greater(momentum / Q[Shortcuts::h], 0.0)) {
19 const auto coulombFriction{mu * cF * Q[Shortcuts::h]};
20 const auto turbulentDrag{momentumSQR * invXi / (Q[Shortcuts::h] * Q[Shortcuts::h])};
21 const auto frictionTerm{-g * (coulombFriction + turbulentDrag) / momentum};
22 S[Shortcuts::hu] = Q[Shortcuts::hu] * frictionTerm;
23 S[Shortcuts::hv] = Q[Shortcuts::hv] * frictionTerm;
24 }
25"""