SPECFEM / specfem2d

SPECFEM2D simulates forward and adjoint seismic wave propagation in two-dimensional acoustic, (an)elastic, poroelastic or coupled acoustic-(an)elastic-poroelastic media, with Convolution PML absorbing conditions.
GNU General Public License v3.0
204 stars 147 forks source link

The Newmark time scheme is poorly accurate in particular in low-velocity fluids #763

Open komatits opened 7 years ago

komatits commented 7 years ago

The Newmark time scheme is poorly accurate in particular in low-velocity fluids. If you need accuracy for that, for now you should likely switch to LDDRK4 (or use Newmark with a small time step), and one day someone should develop a better time scheme (maybe ADER4 as used e.g. by Bruno Lombard.

Analyzed and confirmed by @bottero and @pcristini ; details (in French) below

komatits commented 7 years ago

From @bottero : Le problème c'est le schéma Newmark qui provoque TRES VITE de la dispersion numérique. Essaies en LDDRK tu vas voir.

komatits commented 7 years ago

From @bottero: Oui... ce qui compte c'est la vitesse du milieu. C'est exactement la réflexion que je me suis faite. Je viens de faire d'autres tests... ça se clarifie! J'ai fait les mêmes simus avec du fluide-solide mais avec des vitesses dans le fluide et le solide comparables et j'ai les mêmes problèmes que plus haut : il faut que je passe en LDDRK pour avoir une convergence acceptable ou que je prenne une CFL plus basse (0.2). Donc c'est cohérent! C'est pas un problème de fluide ou de solide (hormis quand on regarde le déplacement parce-qu'on perd un ordre).

Dans le cas ou le modèle est composé de plusieurs milieux avec un des milieux aux fortes disparités de vitesse (comme une couche d'eau et un fond solide dur) on doit baisser le dt si on veut garder la même CFL (afin de s'adapter au milieu qui est plus rapide) et donc on travaille avec un dt très faible pour les calculs dans le milieu lent et ça rattrape l'effet. Par contre si tous les milieux sont rapides apparemment le problème apparaît bien plus tard (comme le benchmark specfem3d globe).

Ma main a couper que si vous aviez fait les mêmes benchmarks pour une terre pleine d'eau vous auriez dû prendre une CFL << 0.3 pour fiter les solutions analytiques avec un schéma d'ordre 2.

Ce qu'il faut retenir c'est qu'il ne suffit pas de régler la CFL pour obtenir la stabilité mais que cela joue FORTEMENT sur la dispersion numérique quand les vitesses sont faibles (c'est le fortement qu'il faut retenir). _Petites vitesses (~vp < 2500m/s) et nb de pas de temps > 5000 => petite CFL (d'autant plus petite que le nb de pas de temps est grand) ou LDDRK _Milieu mixte grandes/petites vitesses => L'effet apparaît plus tard car les zones à haute vitesse font qu'il faut diminuer le DT pour garder une CFL constante. _Milieu avec grandes vitesses (~vp > 2500m/s) et nb de pas de temps < 50000 => pas de problème _Milieu avec grandes vitesses (~vp > 2500m/s) et nb de pas de temps > 50000 => petite CFL (d'autant plus petite que le nb de pas de temps est grand) ou LDDRK

Pour conclure il faut utiliser LDDRK quasi tout le temps quand on simule de l'eau pure sur plusieurs milliers de pas de temps, dans ce cas avec Newmark choisir une CFL > 0.3 est très très mauvais.

komatits commented 7 years ago

Newmark n'est pas super en fait, c'est juste des différences finies centrées d'ordre 2; mais ça marche bien dans les codes FD de Vadim et dans les miens et dans l'industrie pétrolière. Donc y'a pas que l'aspect temps car sinon ça marcherait pas dans les codes FD. Ce qui n'est pas super en fait c'est le fait de mixer de l'ordre élevé spatial avec de l'ordre faible en temps. C'est montré dans De Basabe 2007 je crois de mémoire.

komatits commented 7 years ago

Le livre qui discute ça aussi c'est Hughes 1987 chap 9 mais je sais plus si y'a l'ordre spatial élevé dedans, je crois que non ; mais dans De Basabe 2007 oui.

komatits commented 6 years ago

Just a remark, not really an issue.

ReneSteinmann commented 5 years ago

Hello, I think I might have a similiar problem. I want to model Rayleigh and Love waves in saturated media, which means my poisson's ratio is about 0.49. The sh-simulation works fine, but the p-sv-simulation causes troubles and I have numerical dispersion, probably due to the high difference between s-wave-velocity and p-wave-velocity. Has anyone encountered similar problems or can tell me what is wrong? I also tried LDDRK4 and the newmark time scheme with dt=0.01. Both didn't work.

komatits commented 5 years ago

Hello,

Yes, for media with values of Poisson's ratio very close to 0.5 (say above 0.49 or so) all such classical finite element or spectral element techniques are unsuitable. This is unrelated to the SPECFEM package itself. It is a classical problem in medical imaging for instance, in which P wave speeds are those of water (1500 m/s) and shear wave speeds can be as low as 1 to 2 m/s, leading to Poisson = 0.49999999 . There is no way such very short wave speeds can be discretized accurately at reasonable cost (even in 2D).

See e.g. https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5276 and references therein for a detailed study on two classical ways of addressing that in mechanics and/or medical imaging.

Best wishes, Dimitri.

On 11/14/18 2:30 PM, shirleysheep wrote:

Hello, I think I might have a similar problem. I want to model Rayleigh and Love waves in saturated media, which means my poisson's ratio is about 0.49. The sh-simulation works fine, but the p-sv-simulation causes troubles and I have numerical dispersion, probably due to the high difference between s-wave-velocity and p-wave-velocity. Has anyone encountered similar problems or can tell me what is wrong?

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/geodynamics/specfem2d/issues/763#issuecomment-438662189, or mute the thread https://github.com/notifications/unsubscribe-auth/AFjDKU8VEUzTH0EUC6tks6p-4zXKWO9Iks5uvBsNgaJpZM4NZajE.

-- Dimitri Komatitsch, CNRS Research Director (DR CNRS) Laboratory of Mechanics and Acoustics, Marseille, France http://komatitsch.free.fr

komatits commented 5 years ago

Hello again,

the reference I mentioned below only addresses the first option. The second is e.g. at https://fex.insa-lyon.fr/get?k=HgvxSHZnYI5tbV4Br17 page 147 (link valid for two weeks).

Best wishes, Dimitri.

On 11/21/18 4:39 PM, Dimitri Komatitsch wrote:

Hello,

Yes, for media with values of Poisson's ratio very close to 0.5 (say above 0.49 or so) all such classical finite element or spectral element techniques are unsuitable. This is unrelated to the SPECFEM package itself. It is a classical problem in medical imaging for instance, in which P wave speeds are those of water (1500 m/s) and shear wave speeds can be as low as 1 to 2 m/s, leading to Poisson = 0.49999999 . There is no way such very short wave speeds can be discretized accurately at reasonable cost (even in 2D).

See e.g. https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5276 and references therein for a detailed study on two classical ways of addressing that in mechanics and/or medical imaging.

Best wishes, Dimitri.

On 11/14/18 2:30 PM, shirleysheep wrote:

Hello, I think I might have a similar problem. I want to model Rayleigh and Love waves in saturated media, which means my poisson's ratio is about 0.49. The sh-simulation works fine, but the p-sv-simulation causes troubles and I have numerical dispersion, probably due to the high difference between s-wave-velocity and p-wave-velocity. Has anyone encountered similar problems or can tell me what is wrong?

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/geodynamics/specfem2d/issues/763#issuecomment-438662189, or mute the thread https://github.com/notifications/unsubscribe-auth/AFjDKU8VEUzTH0EUC6tks6p-4zXKWO9Iks5uvBsNgaJpZM4NZajE.

-- Dimitri Komatitsch, CNRS Research Director (DR CNRS) Laboratory of Mechanics and Acoustics, Marseille, France http://komatitsch.free.fr