Adaptive mesh refinement method : numerical density of entropy - - PowerPoint PPT Presentation
Adaptive mesh refinement method : numerical density of entropy production and automatic thresholding. Ersoy, M. 1 a Pons, K. 1 , 2 b 1 EA 2134, IMATH, Universit e de Toulon 2 Principia S.A.S.,13705 La Ciotat cedex 2017, 17 February, LAMA
Adaptive mesh refinement method : numerical density of entropy production and automatic thresholding.
Ersoy, M.1 a Pons, K.1,2 b
1EA 2134, IMATH, Universit´
e de Toulon
2Principia S.A.S.,13705 La Ciotat cedex
2017, 17 February, LAMA – Chamb´ ery
b.
pons@univ-tln.fr
Motivations
numerical simulations of real-life large scale fluid flows in dimension d = 1, 2, 3
(a) Monster waves (Nazare) (b) Wave-breaking (Nazare) (c) Tsunami (Japan) (d) tsunami (Brazil)
AMR 2017, 17 February, LAMA – Chamb´ ery 2 / 22
Motivations
numerical simulations of real-life large scale fluid flows in dimension d = 1, 2, 3 Simulating fluids in large-scale → high memory and computer requirements.
◮ how to compute fast and accurate : → meshing or moving the mesh only in
” desired”regions with a suitable mesh refinement criterion
AMR 2017, 17 February, LAMA – Chamb´ ery 2 / 22
Motivations
numerical simulations of real-life large scale fluid flows in dimension d = 1, 2, 3 Simulating fluids in large-scale → high memory and computer requirements. Mathematical motivations : introducing new tool for numerical purpose !
AMR 2017, 17 February, LAMA – Chamb´ ery 2 / 22
Outline of the talk
Outline of the talk
1 Adaptive Mesh Refinement method
Generality An example
2 Automatic mesh refinement threshold 3 Numerical simulations
AMR 2017, 17 February, LAMA – Chamb´ ery 2 / 22
Outline
Outline
1 Adaptive Mesh Refinement method
Generality An example
2 Automatic mesh refinement threshold 3 Numerical simulations
AMR 2017, 17 February, LAMA – Chamb´ ery 2 / 22
Outline
Outline
1 Adaptive Mesh Refinement method
Generality An example
2 Automatic mesh refinement threshold 3 Numerical simulations
AMR 2017, 17 February, LAMA – Chamb´ ery 3 / 22
Hyperbolic equations and entropy inequality
We focus on general non linear hyperbolic conservation laws ∂w ∂t + ∂f(w) ∂x = 0, (t, x) ∈ R+ × R w(0, x) = w0(x), x ∈ R w ∈ Rd : vector state, f : flux governing the physical description of the flow.
AMR 2017, 17 February, LAMA – Chamb´ ery 4 / 22
Hyperbolic equations and entropy inequality
We focus on general non linear hyperbolic conservation laws ∂w ∂t + ∂f(w) ∂x = 0, (t, x) ∈ R+ × R w(0, x) = w0(x), x ∈ R Weak solutions satisfy S = ∂s(w) ∂t + ∂ψ(w) ∂x = for smooth solution = across rarefaction < across shock where (s, ψ) stands for a convex entropy-entropy flux pair : (∇ψ(w))T = (∇s(w))T Dwf(w)
AMR 2017, 17 February, LAMA – Chamb´ ery 4 / 22
Hyperbolic equations and entropy inequality
We focus on general non linear hyperbolic conservation laws ∂w ∂t + ∂f(w) ∂x = 0, (t, x) ∈ R+ × R w(0, x) = w0(x), x ∈ R Weak solutions satisfy S = ∂s(w) ∂t + ∂ψ(w) ∂x = for smooth solution = across rarefaction < across shock where (s, ψ) stands for a convex entropy-entropy flux pair : (∇ψ(w))T = (∇s(w))T Dwf(w) Entropy inequality ≃“smoothness indicator”
Croisille J.-P., Contribution ` a l’´ Etude Th´ eorique et ` a l’Approximation par ´ El´ ements Finis du Syst` eme Hyperbolique de la Dynamique des Gaz Multidimensionnelle et Multiesp` eces, PhD thesis, Universit´ e de Paris VI, 1991
AMR 2017, 17 February, LAMA – Chamb´ ery 4 / 22
Hyperbolic equations and entropy inequality
We focus on general non linear hyperbolic conservation laws ∂w ∂t + div(f(w)) = 0, (t, x) ∈ R+ × Rd w(0, x) = w0(x), x ∈ Rd Weak solutions satisfy S = ∂s(w) ∂t + div(ψ(w)) = for smooth solution = across rarefaction < across shock where (s, ψ) stands for a convex entropy-entropy flux pair : (∇ψi(w))T = (∇s(w))T Dwfi(w), i = 1, . . . , d Entropy inequality ≃“smoothness indicator”
Croisille J.-P., Contribution ` a l’´ Etude Th´ eorique et ` a l’Approximation par ´ El´ ements Finis du Syst` eme Hyperbolique de la Dynamique des Gaz Multidimensionnelle et Multiesp` eces, PhD thesis, Universit´ e de Paris VI, 1991
AMR 2017, 17 February, LAMA – Chamb´ ery 4 / 22
Finite volume approximation Figure: a cell Ck
Finite volume approximation : wn+1
k
= wn
k − δtn
hk
k+1/2 − F n k−1/2
wn
k ≃ 1
hk
w (tn, x) dx and F n
k+1/2 ≈ 1
δt
f(t, w(t, xk+1/2)) dx
AMR 2017, 17 February, LAMA – Chamb´ ery 5 / 22
Finite volume approximation Figure: a cell Ck
Finite volume approximation : wn+1
k
= wn
k − δtn
hk
k+1/2 − F n k−1/2
wn
k ≃ 1
hk
w (tn, x) dx and F n
k+1/2 ≈ 1
δt
f(t, w(t, xk+1/2)) dx The numerical density of entropy production : Sn
k = sn+1 k
− sn
k
δtn + ψn
k+1/2 − ψn k−1/2
hk
AMR 2017, 17 February, LAMA – Chamb´ ery 5 / 22
Finite volume approximation Figure: a cell Ck
Finite volume approximation : wn+1
k
= wn
k − δtn
hk
F (wn
k, wn a; nk/a)
hk = |Ck|
with wn
k ≃ 1
hk
w (tn, x) dx, and F (wn
k, wn a; nk/a) ≈ 1
δt
f(t, w) · nk/a ds
AMR 2017, 17 February, LAMA – Chamb´ ery 5 / 22
Finite volume approximation Figure: a cell Ck
Finite volume approximation : wn+1
k
= wn
k − δtn
hk
F (wn
k, wn a; nk/a)
hk = |Ck|
with wn
k ≃ 1
hk
w (tn, x) dx, and F (wn
k, wn a; nk/a) ≈ 1
δt
f(t, w) · nk/a ds The numerical density of entropy production : Sn
k = sn+1 k
− sn
k
δtn +
k, wn a; nk/a)
hk
AMR 2017, 17 February, LAMA – Chamb´ ery 5 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k x
2 4 6 8 10
w(x)
1 1.1 1.2 1.3 1.4 1.5
x
2 4 6 8
S(x) = |w′(x)|
0.002 0.004 0.006 0.008 0.01 0.012 S(x) α ¯ S with α = 2
x
2 4 6 8
Refinement indicator
0.2 0.4 0.6 0.8 1 Refined area
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
x
2 4 6 8 10
w(x)
1 1.1 1.2 1.3 1.4 1.5
x
2 4 6 8
S(x) = |w′(x)|
0.002 0.004 0.006 0.008 0.01 0.012 S(x) α ¯ S with α = 2
x
2 4 6 8 0.2 0.4 0.6 0.8 1 Refined area Coarsened area
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ Dynamic mesh refinement : ⋆ Dyadic tree (1D) ⋆ hierarchical numbering : basis 2
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ Dynamic mesh refinement : ⋆ Non-structured grid : macro-cell ⋆ Dyadic tree (1D), Quadtree (2D) ⋆ hierarchical numbering : basis 2,4
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ Dynamic mesh refinement : ⋆ Non-structured grid : macro-cell ⋆ Dyadic tree (1D), Quadtree (2D), Octree (3D) ⋆ hierarchical numbering : basis 2,4,8
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ A simple projection method but the scheme is locally non consistent
[S088, TW05]
Shu C. W., Osher S., Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys., 77(2) :439–471, 1988. Tang H., Warnecke G., A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ A simple projection method but the scheme is locally non consistent
[S088, TW05]
⋆ less time time consuming than other methods Shu C. W., Osher S., Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys., 77(2) :439–471, 1988. Tang H., Warnecke G., A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ A simple projection method but the scheme is locally non consistent
[S088, TW05]
⋆ less time time consuming than other methods ⋆ non consistency is almost negligible if the level of adjacent cells is limited by 2 Shu C. W., Osher S., Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys., 77(2) :439–471, 1988. Tang H., Warnecke G., A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space
Central European Journal of Mathematics, 11(8), pp 1392–1415, 2013.
for three-dimensional two-fluid flows. International Journal of Computational Fluid Dynamics, Taylor & Francis, 2015.
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ A simple projection method but the scheme is locally non consistent
[S088, TW05]
⋆ less time time consuming than other methods ⋆ non consistency is almost negligible if the level of adjacent cells is limited by 2 ◮ Improvement (cpu-time) : local time stepping [EGY13] See more details Shu C. W., Osher S., Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys., 77(2) :439–471, 1988. Tang H., Warnecke G., A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space
Central European Journal of Mathematics, 11(8), pp 1392–1415, 2013.
for three-dimensional two-fluid flows. International Journal of Computational Fluid Dynamics, Taylor & Francis, 2015.
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Mesh refinement indicator : principle & illustration
Compute wn
k
Compute Sn
k : Sn k = 0 =
⇒ the cell is refined or coarsened More precisely : for a given mesh refinement threshold α
◮ Sn
k αS =
⇒ the cell is refined with, for instance, S = 1 |Ω|
Sn
k
◮ Sn
k < αS =
⇒ the cell is coarsened
◮ A simple projection method but the scheme is locally non consistent
[S088, TW05]
⋆ less time time consuming than other methods ⋆ non consistency is almost negligible if the level of adjacent cells is limited by 2 ◮ Improvement (cpu-time) : local time stepping [EGY13] See more details ◮ 2D-3D computer management : BB-AMR [GEYS] See more details Shu C. W., Osher S., Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys., 77(2) :439–471, 1988. Tang H., Warnecke G., A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space
Central European Journal of Mathematics, 11(8), pp 1392–1415, 2013.
for three-dimensional two-fluid flows. International Journal of Computational Fluid Dynamics, Taylor & Francis, 2015.
AMR 2017, 17 February, LAMA – Chamb´ ery 6 / 22
Outline
Outline
1 Adaptive Mesh Refinement method
Generality An example
2 Automatic mesh refinement threshold 3 Numerical simulations
AMR 2017, 17 February, LAMA – Chamb´ ery 7 / 22
An example : the one-dimensional gas dynamics equations for ideal gas
∂ρ ∂t + ∂ρu ∂x = 0 ∂ρu ∂t + ∂
= 0 ∂ρE ∂t + ∂ (ρE + p) u ∂x = 0 p = (γ − 1)ρε where ρ(t, x) : density u(t, x) : velocity p(t, x) : pressure γ := 1.4 : ratio of the specific heats E(ε, u) : total energy ε : internal specific energy E = ε + u2 2
Intel(R) Core(TM) i5-2500 CPU @ 3.30GHz
AMR 2017, 17 February, LAMA – Chamb´ ery 8 / 22
An example : the one-dimensional gas dynamics equations for ideal gas
∂ρ ∂t + ∂ρu ∂x = 0 ∂ρu ∂t + ∂
= 0 ∂ρE ∂t + ∂ (ρE + p) u ∂x = 0 p = (γ − 1)ρε where ρ(t, x) : density u(t, x) : velocity p(t, x) : pressure γ := 1.4 : ratio of the specific heats E(ε, u) : total energy ε : internal specific energy E = ε + u2 2 Conservative variables w = (ρ, ρu, ρE)t convex continuous entropy s(w) = −ρ ln p ργ
Intel(R) Core(TM) i5-2500 CPU @ 3.30GHz
AMR 2017, 17 February, LAMA – Chamb´ ery 8 / 22
Sod’s shock tube problem
Mesh coarsening parameter α : 0.001 , Mesh refinement parameter ¯ S : 1 |Ω|
Sn
kb
CFL : 0.25, Simulation time (s) : 0.4, Initial number of cells : 200, Maximum level of mesh refinement : Lmax .
AMR 2017, 17 February, LAMA – Chamb´ ery 9 / 22
Accuracy
0.2 0.4 0.6 0.8 1 1.2
0.5 1-0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 Density Numerical density of entropy production x ρ on adaptive mesh with Lmax = 4 ρ on uniform fixed mesh N = 681 ρex Sk
n(a) Density and numerical density of en- tropy production.
1 1.5 2 2.5 3 3.5 4 4.5 5
0.5 1 0 0.02 0.04 0.06 0.08 0.1 Mesh refinement level Numerical density of entropy production x level Sk
n|ρ-ρex|
(b) Mesh refinement level, numerical density of entropy production and local error.
Figure: Sod’s shock tube problem : solution at time t = 0.4 s using the AB1M scheme
cells.
AMR 2017, 17 February, LAMA – Chamb´ ery 10 / 22
Properties of S : shock criterion type
Theorem ([P04])
Consider a pth convergent scheme. Let Sn
k be the corresponding numerical
density of entropy production and ∆t = λh be a fixed time step where h stands for the meshsize. Then lim
n→∞ Sn k =
O(∆tp) if the solution is smooth, O 1 ∆t
Puppo G., Numerical entropy production for central schemes. SIAM J. Sci. Comput., 25(4) :1382–1415, 2004.
AMR 2017, 17 February, LAMA – Chamb´ ery 11 / 22
Properties of S : shock criterion type
Theorem ([P04])
Consider a pth convergent scheme. Let Sn
k be the corresponding numerical
density of entropy production and ∆t = λh be a fixed time step where h stands for the meshsize. Then lim
n→∞ Sn k =
O(∆tp) if the solution is smooth, O 1 ∆t
Properties
Consider a monotone scheme. Then, for almost every k, every n, Sn
k 0.
AMR 2017, 17 February, LAMA – Chamb´ ery 11 / 22
Properties of S : shock criterion type
Theorem ([P04])
Consider a pth convergent scheme. Let Sn
k be the corresponding numerical
density of entropy production and ∆t = λh be a fixed time step where h stands for the meshsize. Then lim
n→∞ Sn k =
O(∆tp) if the solution is smooth, O 1 ∆t
Properties
Consider a monotone scheme. Then, for almost every k, every n, Sn
k 0.
Thus, even if locally Sn
k can take positive value, one has Sn k C∆tq,
q p .
See example
AMR 2017, 17 February, LAMA – Chamb´ ery 11 / 22
Outline
Outline
1 Adaptive Mesh Refinement method
Generality An example
2 Automatic mesh refinement threshold 3 Numerical simulations
AMR 2017, 17 February, LAMA – Chamb´ ery 11 / 22
Mesh refinement threshold
Mesh refinement criterion Sn
k and the mesh refinement threshold α
. . . Sn
k αS =
⇒ the cell is refined with S = 1 |Ω|
Sn
k
Sn
k < αS =
⇒ the cell is coarsened . . .
AMR 2017, 17 February, LAMA – Chamb´ ery 12 / 22
Mesh refinement threshold
Mesh refinement criterion Sn
k and the mesh refinement threshold α
. . . Sn
k αS =
⇒ the cell is refined with S = 1 |Ω|
Sn
k
Sn
k < αS =
⇒ the cell is coarsened . . . how to set α suitably ?
AMR 2017, 17 February, LAMA – Chamb´ ery 12 / 22
Mesh refinement threshold
Mesh refinement criterion Sn
k and the mesh refinement threshold α
. . . Sn
k αS =
⇒ the cell is refined with S = 1 |Ω|
Sn
k
Sn
k < αS =
⇒ the cell is coarsened . . . how to set α suitably ? ideally
◮ automatically ◮ accuracy and computational time are“well-balanced” ⋆ catch the smallest maxima but not the smallest
AMR 2017, 17 February, LAMA – Chamb´ ery 12 / 22
Mesh refinement threshold
Mesh refinement criterion Sn
k and the mesh refinement threshold α
. . . Sn
k αS =
⇒ the cell is refined with S = 1 |Ω|
Sn
k
Sn
k < αS =
⇒ the cell is coarsened . . . how to set α suitably ? ideally
◮ automatically ◮ accuracy and computational time are“well-balanced” ⋆ catch the smallest maxima but not the smallest ⋆ simple method (”
low-cost” )
AMR 2017, 17 February, LAMA – Chamb´ ery 12 / 22
Review of classical methods
Mean method or deviation S(x, t) > α 1 |Ω|
S(x, t) dx where α is a tunable dimensionless threshold parameter or deviation method S(x, t) > α 1 |Ω|
S(x, t) dx + βσ(x, t) where σ is the standard deviation with β is a tunable dimensionless threshold parameter
Kallinderis, Y.G., Baron, J.R. : Adaptation methods for a new navier-stokes algorithm. AIAA journal 27(1), 37–43 (1989) Ersoy, M., Golay, F., Yushchenko, L. : Adaptive multiscale scheme based on numerical density of entropy production for conservation laws. Cent.
AMR 2017, 17 February, LAMA – Chamb´ ery 13 / 22
Review of classical methods
Mean method or deviation → pb : set α or (α, β) Filtering : two-steps
x
0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 S(x) α ¯ S with α = 1
(a) Filtering (mean me- thod) : first step
x
0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 min( ¯ S, S(x)) αmin( ¯ S, S(x)) with α = 1
(b) Filtering (mean me- thod) : second step
Aftosmis, M. : Upwind method for simulation of viscous flow on adaptively refined meshes. AIAA journal 32(2), 268–277 (1994) Warren, G.P., Anderson,W.K., Thomas, J.L., Krist, S.L. : Grid convergence for adaptive methods. AIAA paper 1592, 1991 (1991)
AMR 2017, 17 February, LAMA – Chamb´ ery 13 / 22
Review of classical methods
Mean method or deviation → pb : set α or (α, β) Filtering : two-steps → pb : two-step ? Is it enough ? still to set α Filtering : wavelets
x
0.2 0.4 0.6 0.8 1
S(x)
5 10 15 20
(c) S
x
0.2 0.4 0.6 0.8 1
scales
50 100 150
(d) Wavelet coefficients (x = space, y = scale)
Figure: Illustration of the wavelet transformation for a given mesh refinement criterion computed with the Daubechies wavelet with four
vanishing moments (warm colors correspond to large coefficients and cold colors to small coefficients) Leonard, S., Terracol, M., Sagaut, P. : A wavelet-based adaptive mesh refinement criterion for large-eddy simulation. Journal of Turbulence 7(64) (2006)
AMR 2017, 17 February, LAMA – Chamb´ ery 13 / 22
Review of classical methods
Mean method or deviation → pb : set α or (α, β) Filtering : two-steps → pb : two-step ? Is it enough ? still to set α Filtering : wavelets→ pb : efficient but the cost ! still to set α Local maxima→ pb : efficient but the cost ! still to set α
0.0002 0.0004 0.0006 0.0008 0.001 5 10 15 20 25 S(x) x (m)Criterion
(a) S
0.2 0.4 0.6 0.8 1 5 10 15 20 25 Flagged zone(x) x (m)Flagged zone
(b) Region to be refined
Figure: Illustration of the local maxima method for a mesh refinement criterion involving multiple scales
AMR 2017, 17 February, LAMA – Chamb´ ery 13 / 22
Distribution function (L-meas {S(x) > α})
Assumptions and notations
◮ S is smooth and has p local maxima. ◮ S(0) = S(L) = S′(0) = S′(L) = 0 ◮ 0 < S∞ = max
x∈(0,L) S(x) < ∞
x 2 4 6 8 10 S(x) 2 4 6 8 10 12 14 16 18 20 L p = 3 S(x)AMR 2017, 17 February, LAMA – Chamb´ ery 14 / 22
Distribution function (L-meas {S(x) > α})
Assumptions and notations
◮ S is smooth and has p local maxima. ◮ S(0) = S(L) = S′(0) = S′(L) = 0 ◮ 0 < S∞ = max
x∈(0,L) S(x) < ∞
Then
◮
Zα = {x ∈ (0, L); ϕα(x) = S(x) − α = 0 and S′(x) = 0} = ∅ = {x0(α) < x1(α) < · · · < x2pα−2(α) < x2pα−1(α)}
◮ #Zα = 2pα
x 2 4 6 8 10 S(x) 2 4 6 8 10 12 14 16 18 20 L x0(α) x1(α) x2(α) x3(α) pα = 2 p = 3 S(x) αAMR 2017, 17 February, LAMA – Chamb´ ery 14 / 22
Distribution function (L-meas {S(x) > α})
Assumptions and notations
◮ S is smooth and has p local maxima. ◮ S(0) = S(L) = S′(0) = S′(L) = 0 ◮ 0 < S∞ = max
x∈(0,L) S(x) < ∞
Then
◮
Zα = {x ∈ (0, L); ϕα(x) = S(x) − α = 0 and S′(x) = 0} = ∅ = {x0(α) < x1(α) < · · · < x2pα−2(α) < x2pα−1(α)}
◮ #Zα = 2pα ◮ there exist sequence (x∗
k)1kp
x 2 4 6 8 10 S(x) 2 4 6 8 10 12 14 16 18 20 L x0(α) x1(α) x2(α) x3(α) pα = 2 p = 3 x∗ 1 x∗ 2 x∗ 3 S(x) αAMR 2017, 17 February, LAMA – Chamb´ ery 14 / 22
Distribution function (L-meas {S(x) > α})
Assumptions and notations
◮ S is smooth and has p local maxima. ◮ S(0) = S(L) = S′(0) = S′(L) = 0 ◮ 0 < S∞ = max
x∈(0,L) S(x) < ∞
Then
◮
Zα = {x ∈ (0, L); ϕα(x) = S(x) − α = 0 and S′(x) = 0} = ∅ = {x0(α) < x1(α) < · · · < x2pα−2(α) < x2pα−1(α)}
◮ #Zα = 2pα ◮ there exist sequences (x∗
k)1kp and (α∗ k)1kp such that
⋆ ∀k = 1, . . . , p ⋆ S(x∗
k) = α∗ k and S′(x∗ k) = 0, S′′(x∗ k) < 0
⋆ x∗
k /
∈ Zα and α∗
p = S∞
x 2 4 6 8 10 S(x) 2 4 6 8 10 12 14 16 18 20 L p = 3 x∗ 1 x∗ 2 x∗ 3 S(x) α∗ 1 α∗ 2 α∗ 3 = S∞AMR 2017, 17 February, LAMA – Chamb´ ery 14 / 22
Distribution function (L-meas {S(x) > α})
Assumptions and notations
◮ S is smooth and has p local maxima. ◮ S(0) = S(L) = S′(0) = S′(L) = 0 ◮ 0 < S∞ = max
x∈(0,L) S(x) < ∞
With these settings, we define the distribution d α ∈ [0, S∞] → d(α) := 1 if α = 0 , 1 L
pα
x2k+1(α) − x2k(α) if 0 < α < S∞ , if α = S∞ .
AMR 2017, 17 February, LAMA – Chamb´ ery 14 / 22
Properties of d and αP E threshold
Properties
1
S∞ d(α) dα = L S(x) dx.
2 d ∈ C0([0, S∞], R+) and d′ satisfies 1
∀α ∈ [0, S∞], d′(α) < 0
2
∀k ∈ 0, p, lim
α→α∗
k
d′(α) = −∞ with the convention α∗
0 := 0
3 d ∈ Cl
D∗, R+
p−1
(α∗
k, α∗ k+1) .
x 2 4 6 8 10 S(x) 2 4 6 8 10 12 14 16 18 20 L p = 3 x∗ 1 x∗ 2 x∗ 3 S(x) α∗ 1 α∗ 2 α∗ 3 = S∞ 5 10 15 20 d(α) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 α α∗ α∗ 1 α∗ 2 α∗ 3 = S∞ p = 3 d(α)α
5 10 15 20
d′(α)
AMR 2017, 17 February, LAMA – Chamb´ ery 15 / 22
Automatic thresholding : distribution function (L-meas {S(x) > α})
Hard to define α using only d
x
2 4 6 8 10
S(x)
5 10 15 20 25
α
5 10 15 20
d(α)
0.2 0.4 0.6 0.8 1
AMR 2017, 17 February, LAMA – Chamb´ ery 16 / 22
Automatic thresholding : distribution function (L-meas {S(x) > α})
Hard to define α using only d However, one can set α=αD or
◮ either α = αD = min{α | d′(α) = −1}
x
2 4 6 8 10
S(x)
5 10 15 20 25
α
5 10 15 20
d(α)
0.2 0.4 0.6 0.8 1
α
5 10 15 20
d′(α)
Dannenhoffer, J.F. : Grid adaptation for complex two-dimensional transonic flows. Ph.D. thesis, Massachusetts Institute of Technology (1987)
AMR 2017, 17 February, LAMA – Chamb´ ery 16 / 22
Automatic thresholding : distribution function (L-meas {S(x) > α})
Hard to define α using only d However, one can set α=αD or α = αP
◮ either α = αP = min{α | d′′(α) = 0}
x
2 4 6 8 10
S(x)
5 10 15 20 25
α
5 10 15 20
d(α)
0.2 0.4 0.6 0.8 1
α
5 10 15 20
d′(α)
α
5 10 15 20d′′(α)
Powell, K.G., Murman, E.M. : An embedded mesh procedure for leading-edge vortex flows. In : NASA, Langley Research Center, Transonic Symposium : Theory, Application, and Experiment, vol. 1 (1989)
AMR 2017, 17 February, LAMA – Chamb´ ery 16 / 22
Automatic thresholding : distribution function (L-meas {S(x) > α})
Hard to define α using only d However, one can set α=αD or α = αP In practice, d is defined by d(α) =
M−1
dj 1(αj,αj+1)(α) with dj = #{k ; Sn
k > αj} and (αj)0jM =
j M β
0jM
β 1
α
5 10 15 20
d(α)
0.2 0.4 0.6 0.8 1 d(α) uniform grid d(α) non uniform grid
α
5 10 15 20
d′(α)
d′(α) uniform grid filtered d′(α) uniform grid filtered d′(α) non uniform grid Dannenhoffer, J.F. : Grid adaptation for complex two-dimensional transonic flows. Ph.D. thesis, Massachusetts Institute of Technology (1987)
AMR 2017, 17 February, LAMA – Chamb´ ery 16 / 22
Automatic thresholding : distribution function (L-meas {S(x) > α})
Hard to define α using only d However, one can set α=αD or α = αP α is usually too small (too sensitive to β) ! ! ! → detects small fluctuations
AMR 2017, 17 February, LAMA – Chamb´ ery 16 / 22
Automatic thresholding : distribution function (L-meas {S(x) > α})
Hard to define α using only d However, one can set α=αD or α = αP α is usually too small (too sensitive to β) ! ! ! → detects small fluctuations We propose to set α = αP E = min{α ∈ [0, ¯ S] | max{αd(α)}}
α
5 10 15 20
d(α)
0.2 0.4 0.6 0.8 1
α
5 10 15 20
f(α)
0.1 0.2 0.3 0.4 0.5
Pons, K., Ersoy, M. Adaptive mesh refinement method. Part 1 : Automatic thresholding based on a distribution function Pons, K., Ersoy, M. Adaptive mesh refinement method. Part 2 : Application to tsunamis propagation
AMR 2017, 17 February, LAMA – Chamb´ ery 16 / 22
Why αd(α) ?
Properties
Assume that S is twice differentiable and has p local maxima. Then, there exists ∀k = 0 . . . p − 1, α∗∗
k+1 ∈ (α∗ k, α∗ k+1) such that d′′(α∗∗ k+1) = 0
and the function f has p local maxima α1, . . . , αp such that ∀k = 1 . . . p, αk ∈ (α∗∗
k , α∗ k) .
As a consequence αP E > αP > αD.
x 2 4 6 8 10 S(x) 2 4 6 8 10 12 14 16 18 20 L p = 3 x∗ 1 x∗ 2 x∗ 3 S(x) α∗ 1 α∗ 2 α∗ 3 = S∞α
5 10 15 20
d′(α)
AMR 2017, 17 February, LAMA – Chamb´ ery 17 / 22
Illustration : αP E threshold for “discontinuous” flows
α
0.5 1 1.5 2
f(α)
0.02 0.04 0.06 0.08 0.1 0.12
(a) f(α)
x
1 2 3 4 5
S(x)
2 4 6 8 10
S(x) αPE = 0.83974 Sm = 2.4353
(b) S(x)
Figure: The function f for the mesh refinement criterion S(x) = 200 exp(−1000(x − 1.25)2) + 1.25 exp(−5(x − 3.75)2) representing a shock type solution
AMR 2017, 17 February, LAMA – Chamb´ ery 18 / 22
Illustration : αP E threshold for “smooth” flows
α
0.1 0.2 0.3 0.4
f(α)
0.02 0.04 0.06 0.08 0.1 0.12 0.14
(a) f(α)
x
1 2 3 4 5
S(x)
2 4 6 8 10
S(x) αPE = 0.41316 Sm = 0.42152
(b) S(x)
Figure: The function f for the mesh refinement criterion S(x) = 2 exp(−10(x − 1.25)2) + 1.25 exp(−5(x − 3.75)2) representing a smooth flow
AMR 2017, 17 February, LAMA – Chamb´ ery 18 / 22
Outline
Outline
1 Adaptive Mesh Refinement method
Generality An example
2 Automatic mesh refinement threshold 3 Numerical simulations
AMR 2017, 17 February, LAMA – Chamb´ ery 18 / 22
Numerical validation[PE]
A Dam-break problem for the Saint-Venant equations : ∂th + ∂x(hu) = 0, ∂t(hu) + ∂x
where h(t, x) : density u(t, x) : velocity of the water column g : gravity strength Z : topography
K´ evin Pons, Mehmet Ersoy. Adaptive mesh refinement method. Part 1 : Automatic thresholding based on a distribution function
AMR 2017, 17 February, LAMA – Chamb´ ery 19 / 22
Numerical validation[PE]
1 2 3 4 5 6 10 20 30 40 50 60 70 80
0.5 1 1.5 2 2.5 3Water height Mesh refinement criterion x (m) hex(t,x) h(t,x) on adaptive mesh with lmax = 2 Sk
n with α = 0.01265 and Nm = 714Refined zones
(a) A priori error
1 2 3 4 5 6 10 20 30 40 50 60 70 80
50 100 150 200 250 300 350 400Water height Mesh refinement criterion x (m) hex(t,x) h(t,x) on adaptive mesh with lmax = 2 Sk
n with α = 0.05296 and Nm = 667Refined zones
(b) Numerical density of entropy production
1 2 3 4 5 6 10 20 30 40 50 60 70 80
20 40 60 80 100 120 140 160Water height Mesh refinement criterion x (m) hex(t,x) h(t,x) on adaptive mesh with lmax = 2 Sk
n with α = 1.26730 and Nm = 693Refined zones
(c) A posteriori error
1 2 3 4 5 6 10 20 30 40 50 60 70 80
1 2 3 4 5 6 7 8Water height Mesh refinement criterion x (m) hex(t,x) h(t,x) on adaptive mesh with lmax = 2 Sk
n with α = 0.05669 and Nm = 691Refined zones
(d) Gradient of h
Figure: Numerical results for the water height at time t = 2 s
AMR 2017, 17 February, LAMA – Chamb´ ery 19 / 22
Numerical validation[PE]
0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 α(t) t α Sm
(a) A priori error
1 2 3 4 5 6 7 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 α(t) t α Sm
(b) Numerical density of entropy production
0.6 0.8 1 1.2 1.4 1.6 1.8 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 α(t) t α Sm
(c) A posteriori error
0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 α(t) t α Sm
(d) Gradient of h
Figure: Time evolution of the threshold parameter and the mean value Sm
AMR 2017, 17 February, LAMA – Chamb´ ery 19 / 22
Test case I : Solitary wave propagation over a two dimensional reef Figure: Experimental settings and wave gauges locations
K´ evin Pons, Mehmet Ersoy, Fr´ ed´ eric Golay, Richard Marcer. Adaptive mesh refinement method. Part 2 : Application to tsunamis propagation Roeber, V., Cheung, K.F. : Boussinesq-type model for energetic breaking waves in fringing reef environments. Coastal Engineering 70, 1–20 (2012) Roeber, V., Cheung, K.F., Kobayashi, M.H. : Shock-capturing boussinesq-type model for nearshore wave processes. Coastal Engineering 57(4), 407–423 (2010)
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
Adaptive mesh simula- tion Uniform mesh simula- tion Simulation time 240 240 Number of cells 200-560 1000 Re-meshing time step 0.05 s not applicable Time order integration 1 1 Space order integration 1 1 CFL 0.99 0.99
Table: Numerical parameters
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
x[m]
20 40 60 80
time ˜ t
0.2
Exp. h (uniform) h (adaptive) l Z
(a) ˜ t = 55.03
Figure: Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
x[m]
20 40 60 80
time ˜ t
0.2
Exp. h (uniform) h (adaptive) l Z
(a) ˜ t = 66.53
Figure: Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
x[m]
20 40 60 80
time ˜ t
0.2
Exp. h (uniform) h (adaptive) l Z
(a) ˜ t = 69.13
Figure: Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
x[m]
20 40 60 80
time ˜ t
0.2
Exp. h (uniform) h (adaptive) l Z
(a) ˜ t = 70.68
Figure: Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
x[m]
20 40 60 80
time ˜ t
0.2
Exp. h (uniform) h (adaptive) l Z
(a) ˜ t = 76.33
Figure: Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
x[m]
20 40 60 80
time ˜ t
0.2
Exp. h (uniform) h (adaptive) l Z
(a) ˜ t = 125.03
Figure: Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
time ˜ t
100 150 200 η h0
0.1 0.2 0.3 0.4 Experiment Numerical (uniform) Numerical (adaptive)
(a) x = 17.6m
time ˜ t
100 150 200 η h0
0.1 0.2 0.3 0.4 Experiment Numerical (uniform) Numerical (adaptive)
(b) x = 50.4m
time ˜ t
100 150 200 η h0
0.1 0.2 0.3 0.4 Experiment Numerical (uniform) Numerical (adaptive)
(c) x = 58.1m
time ˜ t
100 150 200 η h0
0.1 0.2 0.3 0.4 Experiment Numerical (uniform) Numerical (adaptive)
(d) x = 65.2m
Figure: Surface profiles of solitary wave propagation in time at wave gauges 1 to 6.
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case I : Solitary wave propagation over a two dimensional reef
0.01 0.02 0.03 0.04 0.05 0.06 0.07 60 80 100 120 140 160 180 200 220 α time t ~ Threshold Sm
(a) Threshold
200 400 600 800 1000 1200 60 80 100 120 140 160 180 200 220 Number of cells time t ~ Adaptive mesh Uniform mesh
(b) Number of cells
Figure: Time evolution of the mesh refinement threshold and the number of cells : 95 s against 210 s
AMR 2017, 17 February, LAMA – Chamb´ ery 20 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf Figure: Experimental settings
K´ evin Pons, Mehmet Ersoy, Fr´ ed´ eric Golay, Richard Marcer. Adaptive mesh refinement method. Part 2 : Application to tsunamis propagation. Lynett, P.J., Swigler, D., Son, S., Bryant, D., Socolofsky, S. : Experimental study of solitary wave evolution over a 3d shallow shelf. Coastal Engineering Proceedings 1(32), 1 (2011).
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
Adaptive mesh simula- tion Uniform mesh simula- tion Simulation time 30 s 30 s Number of blocks 128 128 Number of cells 7 500-25 000 33 000 Re-meshing time step 0.25 s not applicable Time order integration 2 2 Space order integration 2 2 CFL 0.5 0.5
Table: Numerical parameters
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
(a) t=0.5s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
(a) t=2.5s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
(a) t=5.75s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
(a) t=23.75s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 5 10 15 20 25 30 35 Free surface [m] time [s] Experiment Numerical (uniform) Numerical (adaptive)
(a) WG2 results
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 5 10 15 20 25 30 35 Free surface [m] time [s] Experiment Numerical (uniform) Numerical (adaptive)
(b) WG4 results
0.05 0.1 0.15 0.2 5 10 15 20 25 30 35 Free surface [m] time [s] Experiment Numerical (uniform) Numerical (adaptive)
(c) WG6 results
0.05 0.1 0.15 0.2 0.25 5 10 15 20 25 30 35 Free surface [m] time [s] Experiment Numerical (uniform) Numerical (adaptive)
(d) WG7 results
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case II : Solitary wave propagation over an irregular 3-d shallow shelf
0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 5 10 15 20 25 30 35 40 45 α time [s] Threshold Sm
(a) Threshold
5000 10000 15000 20000 25000 30000 35000 40000 5 10 15 20 25 30 35 40 45 50 Number of cells time [s] Adaptive mesh Uniform mesh
(b) Number of cells
Figure: Time evolution of the mesh refinement threshold and the number of cells : speed up the computation by 2.5 time
AMR 2017, 17 February, LAMA – Chamb´ ery 21 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
(a) Top view (b) Side view
Figure: Settings
K´ evin Pons, Mehmet Ersoy, Fr´ ed´ eric Golay, Richard Marcer. Adaptive mesh refinement method. Part 2 : Application to tsunamis propagation
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
Adaptive mesh simula- tion Uniform mesh simula- tion Simulation time 30 s 30 s Number of blocks 240 240 Number of cells 8 000-40 000 62 000 Re-meshing time step 0.25 s not applicable Time order integration 2 2 Space order integration 1 1 CFL 0.5 0.5
Table: Numerical parameters
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
(a) t = 11.25 s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
(a) t = 13.25 s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
(a) t = 16 s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
(a) t = 17.5 s
Figure: Numerical water height (coloration is issue from the kinetic energy)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
(a) Gauge 1
(b) Gauge 2
(c) Gauge 3
Figure: Free surface results at different positions : experimental data versus numerical simulation with and without mesh adaptivity
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Test case III : Tsunami runup onto a complex three dimensional Monai-Walley beach
5x10-5 0.0001 0.00015 0.0002 5 10 15 20 25 30 α time [s] Threshold Sm
(a) Threshold
10000 20000 30000 40000 50000 60000 70000 5 10 15 20 25 30 Number of cells time [s] Adaptive mesh Uniform mesh
(b) Number of cells
Figure: Time evolution of the mesh refinement threshold and the number of cells : speed up the computation by 3 time
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Example
Let us consider the transport equation : wt + wx = w(0, x) = w0(x)
Back to slide
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Example
Let us consider the transport equation : wt + wx = w(0, x) = w0(x) and the Godunov scheme : wn+1
k
= wn
k − δt
δx
k − wn k−1
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Example
Let us consider the transport equation : wt + wx = w(0, x) = w0(x) and the Godunov scheme : wn+1
k
= wn
k − δt
δx
k − wn k−1
k
= s(wn+1
k
) − s(wn
k)
δt + ψ(s(wn
k)) − ψ(s(wn k−1))
δx with s(w) = w2 and ψ(w) = w2.
Back to slide
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Example
Let us consider the transport equation : wt + wx = w(0, x) = w0(x) and the Godunov scheme : wn+1
k
= wn
k − δt
δx
k − wn k−1
k
= s(wn+1
k
) − s(wn
k)
δt + ψ(s(wn
k)) − ψ(s(wn k−1))
δx with s(w) = w2 and ψ(w) = w2. Substituting wn+1
k
into Sn+1
k
, we get Sn+1
k
= −ε wn
k − wn k−1
δx 2 0 with ε = δx
δx
Back to slide
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Time restriction
Explicit adaptive schemes : time consuming due to the restriction wδt h 1, h = min
k hk
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Time restriction, local time stepping approach
Explicit adaptive schemes : time consuming due to the restriction wδt h 1, h = min
k hk
Local time stepping algorithm : save the cpu-time
◮ Sort cells in groups w.r.t. to their level Muller S., Stiriba Y., Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping. SIAM J. Sci. Comput., 30(3) :493 ?531, 2007.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Time restriction, local time stepping approach & Aims
Explicit adaptive schemes : time consuming due to the restriction wδt h 1, h = min
k hk
Local time stepping algorithm : save the cpu-time
◮ Sort cells in groups w.r.t. to their level ◮ Update the cells following the local time stepping algorithm. Muller S., Stiriba Y., Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping. SIAM J. Sci. Comput., 30(3) :493 ?531, 2007.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Illustration Figure: t = tn
with δF n
k−1,k,k+1 :=
k+1/2(wk, wk+1) − F n k−1/2(wk−1, wk)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Illustration Figure: tn1 = tn + δtn
wn1
k000 = wn k000 − δtn
hk000 δF n
k00,k000,k001
wn1
k001 = wn k001 − δtn
hk001 δF n
k000,k001,k+1b
with δF n
k−1,k,k+1 :=
k+1/2(wk, wk+1) − F n k−1/2(wk−1, wk)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Illustration Figure: tn2 = tn + 2δtn
wn2
k00 = wn1 k00 − δtn
hk00 δF n1
k−10,k00,k000
wn2
k000 = wn1 k000 − δtn
hk000 δF n1
k00,k000,k001
wn2
k001 = wn1 k001 − δtn
hk001 δF n1
k000,k001,k+1b
with δF n
k−1,k,k+1 :=
k+1/2(wk, wk+1) − F n k−1/2(wk−1, wk)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Illustration Figure: tn3 = tn + 3δtn
wn3
k000 = wn2 k000 − δtn
hk000 δF n2
k00,k000,k001
wn3
k001 = wn2 k001 − δtn
hk001 δF n2
k000,k001,k+1b
with δF n
k−1,k,k+1 :=
k+1/2(wk, wk+1) − F n k−1/2(wk−1, wk)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Illustration Figure: tn+1 = tn + 4δtn
wn+1
k−10 = wn3 k−10 −
δtn hk−10 δF n3
k−2b,k−10,k00
wn+1
k00 = wn3 k00 − δtn
hk00 δF n3
k−10,k00,k000
wn+1
k000 = wn3 k000 − δtn
hk000 δF n3
k00,k000,k001
wn+1
k001 = wn3 k001 − δtn
hk001 δF n3
k000,k001,k+1b
with δF n
k−1,k,k+1 :=
k+1/2(wk, wk+1) − F n k−1/2(wk−1, wk)
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
local time stepping algorithm
foreach i ∈ {1, 2N} do Let j be the biggest integer such that 2j divides i foreach interface xk+1/2 such that Lk+1/2 N − j do
1 compute the integral of Fk+1/2(t) on the time interval 2N−Lk+1/2δtn, 2 distribute Fk+1/2(tn) to the two adjacent cells, 3 update only the cells of level greater than N − j.
end end
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Shu and Osher test case
Mesh coarsening parameter α : 0.001 , Mesh refinement parameter ¯ S : 1 |Ω|
Sn
kb
CFL : 0.219, Simulation time (s) : 0.18, Initial number of cells : 500, Maximum level of mesh refinement : Lmax = 4.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Efficiency of the local time stepping method
P ρ − ρrefl1
x
cpu-time NLmax maximum number of cells AB1 0.288 4.74 10−2 181 1574 2308 AB1M 0.288 4.80 10−2 120 1572 2314
Table: Shu and Osher test case : comparison of numerical schemes of order 1
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
P ρ − ρrefl1
x
cpu-time NLmax maximum number of cells AB1 0.288 4.74 10−2 181 1574 2308 AB1M 0.288 4.80 10−2 120 1572 2314 AB2 0.287 2.75 10−2 170 1391 2023
Table: Shu and Osher test case : comparison of numerical schemes of order 1 and 2
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
P ρ − ρrefl1
x
cpu-time NLmax maximum number of cells AB1 0.288 4.74 10−2 181 1574 2308 AB1M 0.288 4.80 10−2 120 1572 2314 AB2 0.287 2.75 10−2 170 1391 2023 AB2M 0.286 2.74 10−2 108 1357 1994
Table: Shu and Osher test case : comparison of numerical schemes of order 1 and 2
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
P ρ − ρrefl1
x
cpu-time NLmax maximum number of cells AB1 0.288 4.74 10−2 181 1574 2308 AB1M 0.288 4.80 10−2 120 1572 2314 AB2 0.287 2.75 10−2 170 1391 2023 AB2M 0.286 2.74 10−2 108 1357 1994 RK2 0.285 2.08 10−2 299 1375 2005
Table: Shu and Osher test case : comparison of numerical schemes of order 1 and 2
European Journal of Mathematics, 11(8), pp 1392-1415, 2013.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
P ρ − ρrefl1
x
cpu-time NLmax maximum number of cells AB1 0.288 4.74 10−2 181 1574 2308 AB1M 0.288 4.80 10−2 120 1572 2314 AB2 0.287 2.75 10−2 170 1391 2023 AB2M 0.286 2.74 10−2 108 1357 1994 RK2 0.285 2.08 10−2 299 1375 2005
Table: Shu and Osher test case : comparison of numerical schemes of order 1 and 2
European Journal of Mathematics, 11(8), pp 1392-1415, 2013.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Reference solution&Numerical results
(a) Density and numerical density of en- tropy production.
3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 Density x AB1 AB2 RK2 Reference solution
(b) Zoom on oscillating region.
Figure: Shu and Osher test case.
Back to slide
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. For fast computation, the following are required
◮ parallel treatment ◮ hierarchical grids
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. Some interesting issues :
◮ 2D quad-tree [1], ◮ 3D octree [2], ◮ 2/3D anisotropic AMR [3]. Zhang, M., and W.M. Wu. 2011. A two dimensional hydrodynamic and sediment transport model for dam break based on finite volume method with quadtree grid. Applied Ocean Research 33 (4) : 297 – 308. Losasso, F., F. Gibou, and R. Fedkiw. 2004. Simulating Water and Smoke with an Octree Data Structure. ACM Trans. Graph. 23 (3) : 457–462, 2004. Hachem, E., S. Feghali, R. Codina, and T. Coupez. Immersed stress method for fluid structure interaction using anisotropic mesh adaptation. International Journal for Numerical Methods in Engineering 94 (9) : 805–825, 2013.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. The strategy adopted :
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. The strategy adopted :
1
1 domain= 1 block=1 cpu :“failure”→ synchronization depends on the finest domain
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. The strategy adopted :
1
1 domain= 1 block=1 cpu :“failure”→ synchronization depends on the finest domain
2
1 domain= n × blocks = 1cpu :“good compromise”→ each domain has almost the same number number of cells →“better”synchronization
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. The strategy adopted :
1
1 domain= 1 block=1 cpu :“failure”→ synchronization depends on the finest domain
2
1 domain= n × blocks = 1cpu :“good compromise”→ each domain has almost the same number number of cells →“better”synchronization
3
It certainly exists better strategy . . .
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
two and three dimensional case : BB-AMR
Main difficulty : mesh and data structure. The strategy adopted : Management of domain’s interfaces, projection step, . . .
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering more sophisticated numbering exists . . .
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering more sophisticated numbering exists . . . re-numbering and re-meshing being expensive
◮ the mesh should kept constant on a time interval
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering more sophisticated numbering exists . . . re-numbering and re-meshing being expensive
◮ the mesh should kept constant on a time interval ◮ AMR time-step computed through the smallest block and not the smallest cell
Tn+1 − Tn = ∆TAMR is given by the CFL ∆TAMR β mink hblockk maxk ublockk, 0 < β 1.
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
domain= n × blocks = 1cpu
How it works ? each domain has almost the same number of cells domain are defined using Cuthill-McKee numbering more sophisticated numbering exists . . . re-numbering and re-meshing being expensive
◮ the mesh should kept constant on a time interval ◮ AMR time-step computed through the smallest block and not the smallest cell ◮ Gain is important and numerical stability is conserved !
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Examples :
A two dimensional example of BB-AMR with 3 domains and 9 blocks.
(a) AMR T0 (b) AMR T1 (c) AMR T2
Back to slide
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22
Examples :
A two dimensional example of BB-AMR with 3 domains and 9 blocks.
(f) AMR T0 (g) AMR T1 (h) AMR T2
A three dimensional example of BB-AMR with 3 domains and 27 blocks.
(i) Block-based mesh (j) Domain decomposition
Back to slide
AMR 2017, 17 February, LAMA – Chamb´ ery 22 / 22