Improvements to a semi-implicit integration scheme for viscous fluids in SPH
Giuseppe BilottaINGV Vito ZagoNWU Veronica CentorrinoINGV Alexis HéraultCNAM INGV Robert A. DalrympleNWU Ciro Del NegroINGV SPHERIC 2021
INGV Istituto Nazionale di Geofisica e Vulcanologia
Low Reynolds number ⇒ \Longrightarrow
Standard WCSPH with explicit integrator:
Δ t ≤ min β { C 1 h ∥ a ⃗ β ∥ , C 2 h c β , C 3 ρ β h 2 μ β } \Delta t \le \min_\beta \left\{
C_1\sqrt{\frac{h}{\lVert\vec a_\beta\rVert}},
C_2\frac{h}{c_\beta},
C_3\frac{\rho_\beta h^2}{\mu_\beta}
\right\}
h h c c a ⃗ \vec a ρ \rho μ \mu
High viscosity + moderate/high resolution ⇒ \Longrightarrow Δ t \Delta t
(For sound speed we cheat with c = 10 u m a x c = 10 u_{max}
Implicit integration to the rescue
Avoid viscous constraint on time-step with semi -implicit integration (Zago et al. , 2018 doi:10.1016/j.jcp.2018.07.060 ):
u ⃗ β ( n + 1 ) = u ⃗ β ( n ) + Δ t ( f ⃗ β , P ( n ) + g ⃗ ) + Δ t ∑ α 2 μ ‾ α β ( n ) ρ α ( n ) ρ β ( n ) F α β ( n ) m α u ⃗ α β ( n + 1 )
\vec u_\beta^{(n+1)} = \vec u_\beta^{(n)} +
\Delta t \left(\vec f_{\beta,P}^{(n)} + \vec g\right) +
\Delta t \sum_\alpha
\frac{2\bar\mu_{\alpha\beta}^{(n)}}{\rho_\alpha^{(n)} \rho_\beta^{(n)}}
F^{(n)}_{\alpha\beta} m_\alpha \vec u_{\alpha\beta}^{(n+1)}
u u F F ∇ W \nabla W f ⃗ P \vec f_P ⋅ ( ⋅ ) \cdot^{(\cdot)}
Rearrange and get A ( n ) ( Δ t ) u ⃗ ( n + 1 ) = b ⃗ ( n ) ( Δ t )
A^{(n)}(\Delta t) \vec u^{(n+1)} = \vec b^{(n)}(\Delta t)
b ⃗ \vec b
Simple case
Fluid part of matrix A ( n ) A^{(n)}
Boundary? Depends on BC!
Dynamic BC ⇒ \Longrightarrow ⇒ \Longrightarrow ⇒ \Longrightarrow
(But let’s be honest here, there’s better BCs.)
Other boundary conditions?
Consider e.g. “dummy” boundary conditions (Adami et al. ).
Viscous contribution from boundary 𝒲 \mathcal W to fluid ℱ \mathcal F
u ⃗ β , v = 2 u ⃗ β , w − ∑ α ∈ ℱ u ⃗ α W β α ∑ α ∈ ℱ W β α
\vec u_{\beta,v} = 2 \vec u_{\beta,w} - \frac{
\sum_{\alpha\in\mathcal F} \vec u_\alpha W_{\beta\alpha}
}{
\sum_{\alpha\in\mathcal F} W_{\beta\alpha}
}
u ⃗ w \vec u_w
Implicit visc. scheme ⇒ \Longrightarrow u ⃗ v \vec u_v 𝒲 \mathcal W u ⃗ \vec u ℱ \mathcal F together .
Matrix shape
Diagonal:
Off-diagonal
A β β ( n ) ( Δ t ) = { 1 − Δ t ∑ α ∈ ℱ ∪ 𝒲 K α β ( n ) ∀ β ∈ ℱ , 1 ∀ β ∈ 𝒲
A_{\beta\beta}^{(n)}(\Delta t) = \begin{cases}
1 - \Delta t \displaystyle\sum_{\alpha\in\mathcal F \cup \mathcal W} K_{\alpha\beta}^{(n)}
&\forall\beta\in\mathcal F,\\
1 & \forall\beta\in\mathcal W
\end{cases}
A β α ( n ) ( Δ t ) = { Δ t K α β ( n ) β ∈ ℱ , neib α ∈ ℱ ∪ 𝒲 W β α ∑ α ′ ∈ ℱ W β α ′ β ∈ 𝒲 , neib α ∈ ℱ 0 otherwise
A_{\beta\alpha}^{(n)}(\Delta t) = \begin{cases}
\Delta t K_{\alpha\beta}^{(n)} &
\text{$\beta\in\mathcal F$, neib $\alpha \in\mathcal F\cup\mathcal W$}\\
\frac{W_{\beta\alpha}}{\displaystyle\sum_{\alpha'\in\mathcal F} W_{\beta\alpha'}}
& \text{$\beta\in\mathcal W$, neib $\alpha\in\mathcal F$}\\
0 & \text{otherwise}
\end{cases}
with K α β ( n ) = − 2 μ ‾ α β ( n ) ρ α ( n ) ρ β ( n ) F α β ( n ) m α K_{\alpha\beta}^{(n)} = -\frac{2\bar\mu_{\alpha\beta}^{(n)}}{\rho_\alpha^{(n)} \rho_\beta^{(n)}} F_{\alpha\beta}^{(n)} m_\alpha
Fluid rows still diagonally dominant: | A β β | = 1 + ∑ α ≠ β | A β α | \lvert A_{\beta\beta}\rvert = 1 + \sum_{\alpha \ne \beta} \lvert A_{\beta\alpha}\rvert not (| A β β | = ∑ α ≠ β | A β α | \lvert A_{\beta\beta}\rvert = \sum_{\alpha \ne \beta} \lvert A_{\beta\alpha}\rvert
Matrix is Weakly-Chained Diagonally Dominant (chain: single link from 𝒲 \mathcal W ℱ \mathcal F ⇒ \Longrightarrow
Can still solve (?)
Matrix is non-singular, but also not symmetric.A ( n ) ( Δ t ) u ⃗ ( n + 1 ) = b ⃗ ( n ) ( Δ t ) A^{(n)}(\Delta t) \vec u^{(n+1)} = \vec b^{(n)}(\Delta t)
In single precision.
With residue tolerance ε t o l ∥ b ⃗ ∥ \varepsilon_{tol} \lVert \vec b\rVert
And ε t o l = ε M \varepsilon_{tol} = \varepsilon_M
(Narrator voice: they couldn’t)
BiCGSTAB numerical stability
BiCGSTAB not STAB enough?
(Everybody knows that: it’s why they use ε t o l = 10 − 5 \varepsilon_{tol} = 10^{-5}
But why go the easy way when you can suffer?
Let’s fix BiCGSTAB instead.
Standard BiCGSTAB
x ⃗ 0 \vec x_0 r ⃗ 0 = b ⃗ − A x ⃗ 0 \vec r_0 = \vec b - A\vec x_0 r ⃗ ̂ 0 = r ⃗ 0 \hat{\vec r}_0 = \vec r_0 p ⃗ 0 = 0 ⃗ \vec p_0 = \vec 0 γ 0 = α 0 = ω 0 = 1 \gamma_0 = \alpha_0 = \omega_0 = 1
γ i = r ⃗ ̂ 0 ⋅ r ⃗ i − 1 \gamma_i = \hat{\vec r}_0 \cdot \vec r_{i-1} β i = γ i γ i − 1 α i − 1 ω i − 1 \beta_i = \frac{\gamma_i}{\gamma_{i-1}} \frac{\alpha_{i-1}}{\omega_{i-1}} p ⃗ i = r ⃗ i − 1 + β i ( p ⃗ i − 1 − ω i − 1 A p ⃗ i − 1 ) \vec p_i = \vec r_{i-1} + \beta_i\left( \vec p_{i-1} - \omega_{i-1} A \vec p_{i-1} \right) δ i = r ⃗ ̂ 0 ⋅ A p ⃗ i \delta_i = \hat{\vec r}_0 \cdot A \vec p_i α i = γ i δ i \alpha_i = \frac{\gamma_i}{\delta_i} r ⃗ ⋆ = r ⃗ i − 1 − α i A p ⃗ i \vec r_\star = \vec r_{i-1} - \alpha_i A \vec p_i x ⃗ ⋆ = x ⃗ i − 1 + α i p ⃗ i \vec x_\star = \vec x_{i-1} + \alpha_i \vec p_i ω i = r ⃗ ⋆ ⋅ A r ⃗ ⋆ ( A r ⃗ ⋆ ) ⋅ ( A r ⃗ ⋆ ) \omega_i = \frac{\vec r_\star \cdot A \vec r_\star}{(A \vec r_\star)\cdot(A \vec r_\star)} r ⃗ i = r ⃗ ⋆ − ω i A r ⃗ ⋆ \vec r_i = \vec r_\star - \omega_i A \vec r_\star x ⃗ i = x ⃗ ⋆ + ω i r ⃗ ⋆ \vec x_i = \vec x_\star + \omega_i \vec r_\star
(new ⇒)
Improving BiCGSTAB/I
x ⃗ 0 \vec x_0 r ⃗ 0 = b ⃗ − A x ⃗ 0 \vec r_0 = \vec b - A\vec x_0 r ⃗ ̂ 0 = r ⃗ 0 \hat{\vec r}_0 = \vec r_0 p ⃗ 0 = 0 ⃗ \vec p_0 = \vec 0 γ 0 = α 0 = ω 0 = 1 \gamma_0 = \alpha_0 = \omega_0 = 1
γ i = r ⃗ ̂ 0 ⋅ r ⃗ i − 1 \gamma_i = \hat{\vec r}_0 \cdot \vec r_{i-1} β i = γ i γ i − 1 α i − 1 ω i − 1 \beta_i = \frac{\gamma_i}{\gamma_{i-1}} \frac{\alpha_{i-1}}{\omega_{i-1}} p ⃗ i = r ⃗ i − 1 + β i ( p ⃗ i − 1 − ω i − 1 A p ⃗ i − 1 ) \vec p_i = \vec r_{i-1} + \beta_i\left( \vec p_{i-1} - \omega_{i-1} A \vec p_{i-1} \right) ⇐ \Longleftarrow δ i = r ⃗ ̂ 0 ⋅ A p ⃗ i \delta_i = \hat{\vec r}_0 \cdot A \vec p_i α i = γ i δ i \alpha_i = \frac{\gamma_i}{\delta_i} r ⃗ ⋆ = r ⃗ i − 1 − α i A p ⃗ i \vec r_\star = \vec r_{i-1} - \alpha_i A \vec p_i x ⃗ ⋆ = x ⃗ i − 1 + α i p ⃗ i \vec x_\star = \vec x_{i-1} + \alpha_i \vec p_i ω i = r ⃗ ⋆ ⋅ A r ⃗ ⋆ ( A r ⃗ ⋆ ) ⋅ ( A r ⃗ ⋆ ) \omega_i = \frac{\vec r_\star \cdot A \vec r_\star}{(A \vec r_\star)\cdot(A \vec r_\star)} r ⃗ i = r ⃗ ⋆ − ω i A r ⃗ ⋆ \vec r_i = \vec r_\star - \omega_i A \vec r_\star x ⃗ i = x ⃗ ⋆ + ω i r ⃗ ⋆ \vec x_i = \vec x_\star + \omega_i \vec r_\star
Improving BiCGSTAB/II
Expand and simplify! Consider: p ⃗ i = r ⃗ i − 1 + β i ( p ⃗ i − 1 − ω i − 1 A p ⃗ i − 1 ) = r ⃗ i − 1 + β i p ⃗ i − 1 − β i ω i − 1 A p ⃗ i − 1 \vec p_i = \vec r_{i-1} + \beta_i \left( \vec p_{i-1} - \omega_{i-1} A \vec p_{i-1} \right) = \vec r_{i-1} + \beta_i \vec p_{i-1} - \beta_i \omega_{i-1} A \vec p_{i-1} β i = γ i γ i − 1 α i − 1 ω i − 1 = γ i γ i − 1 γ i − 1 δ i − 1 ω i − 1 = γ i δ i − 1 ω i − 1 = α ′ i − 1 ω i − 1
\beta_i =
\frac{\gamma_i}{\gamma_{i-1}} \frac{\alpha_{i-1}}{\omega_{i-1}} =
\frac{\gamma_i}{\gamma_{i-1}} \frac{\gamma_{i-1}}{\delta_{i-1}\omega_{i-1}} =
\frac{\gamma_i}{\delta_{i-1}\omega_{i-1}} =
\frac{\alpha'_{i-1}}{\omega_{i-1}}
α ′ i − 1 = γ i δ i − 1
\alpha'_{i-1} = \frac{\gamma_i}{\delta_{i-1}}
p ⃗ i = r ⃗ i − 1 + β i p ⃗ i − 1 − α ′ i − 1 A p ⃗ i − 1 \vec p_i = \vec r_{i-1} + \beta_i \vec p_{i-1} - \alpha'_{i-1} A \vec p_{i-1}
Modified BiCGSTAB
x ⃗ 0 , r ⃗ 0 = b ⃗ − A ⃗ x ⃗ 0 , r ⃗ ̂ 0 = r ⃗ 0 \vec x_0, \vec r_0 = \vec b - \vec A\vec x_0, \hat{\vec r}_0 = \vec r_0 p ⃗ 0 \vec p_0 γ 0 = r ⃗ ̂ 0 ⋅ r ⃗ 0 , α ′ 0 = β 0 = 0 \gamma_0 = \hat{\vec r}_0 \cdot \vec r_0, \alpha'_0 = \beta_0 = 0
p ⃗ i = r ⃗ i − 1 + β i − 1 p ⃗ i − 1 − α ′ i − 1 A p ⃗ i − 1 \vec p_i = \vec r_{i-1} + \beta_{i-1} \vec p_{i-1} - \alpha'_{i-1} A \vec p_{i-1} δ i = r ⃗ ̂ 0 ⋅ A p ⃗ i \delta_i = \hat{\vec r}_0 \cdot A \vec p_i α i = γ i δ i \alpha_i = \frac{\gamma_i}{\delta_i} r ⃗ ⋆ = r ⃗ i − 1 − α i A p ⃗ i \vec r_\star = \vec r_{i-1} - \alpha_i A \vec p_i x ⃗ ⋆ = x ⃗ i − 1 + α i p ⃗ i \vec x_\star = \vec x_{i-1} + \alpha_i \vec p_i ω i = r ⃗ ⋆ ⋅ A r ⃗ ⋆ ( A r ⃗ ⋆ ) ⋅ ( A r ⃗ ⋆ ) \omega_i = \frac{\vec r_\star \cdot A \vec r_\star}{(A \vec r_\star)\cdot(A \vec r_\star)} r ⃗ i = r ⃗ ⋆ − ω i A r ⃗ ⋆ \vec r_i = \vec r_\star - \omega_i A \vec r_\star x ⃗ i = x ⃗ ⋆ + ω i r ⃗ ⋆ \vec x_i = \vec x_\star + \omega_i \vec r_\star γ i + 1 = r ⃗ ̂ 0 ⋅ r ⃗ i \gamma_{i+1} = \hat{\vec r}_0 \cdot \vec r_i α ′ i = γ i + 1 δ i \alpha'_i = \frac{\gamma_{i+1}}{\delta_i} β i = α ′ i ω i \beta_i = \frac{\alpha'_i}{\omega_i}
(⇐ old)
Non-Newtonian Poiseuille
Poiseuille w/ Bingham rheology (Papanastasiou regularization),
Multi-GPU scaling
Implementation is multi-GPU-capable (easier with reworked BiCGSTAB).
INGVIstituto Nazionale di Geofisica e Vulcanologia
