[cigcommits] r19920  seismo/3D/SPECFEM3D_GLOBE/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Apr 4 15:22:43 PDT 2012
Author: dkomati1
Date: 20120404 15:22:43 0700 (Wed, 04 Apr 2012)
New Revision: 19920
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt
Log:
updated the todo list
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt
===================================================================
 seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt 20120404 16:09:11 UTC (rev 19919)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/todo_list_please_dont_remove.txt 20120404 22:22:43 UTC (rev 19920)
@@ 4,117 +4,33 @@
Things that could be done in a future version:
+ fix the 3D attenuation memory size problem by making it optional and off by default.
+To do that, the following command will be very useful:
+nm printsize sizesort reversesort radix=d ./bin/xspecfem3D  cut f 2 d ' '  less
+
 check if sigma_xz is inverted with sigma_zx and so on in compute_forces; usually this does not matter because the stress tensor is symmetric, but for Cosser
at or for some CPML formulations this could matter
 use a more accurate time scheme: either symplectic schemes (but there could be a problem in the viscoelastic case, see below) or else Runge Kutta RK4, which is non staggered and would thus solve the approximate implementation of the displacement gradient that I used to implement attenuation back in 1999. It is not clear if symplectic schemes could be suitable for viscoelastic media (Tarje implemented them for elastic media only), but RK4 is definitely suitable (Roland Martin used it successfully in 2D viscoelastic).
+ implement RK4 and LDDRK46 highorder time schemes, including for viscoelastic media: will be done by Zhinan Xie
Subject: Re: SEM attenuation
Date: Fri, 23 Jul 2010 12:16:40 +0200
From: Tarje NissenMeyer <tarjen at ethz.ch>
Organization: ETH Zurich
To: Dimitri Komatitsch <dimitri.komatitsch at univpau.fr>
CC: <yingz at vt.edu>, Jeroen Tromp <jtromp at princeton.edu>, Min Chen <mchen at gps.caltech.edu>, Vala Hjorleifsdottir <vala at ldeo.columbia.edu>, "Brian Savage" <savage at uri.edu>, ShiannJong Lee <sjlee at earth.sinica.edu.tw>, "Roland Martin" <roland.martin at univpau.fr>, Bernhard Schuberth <mail at bernhardschuberth.de>, Carl Tape <carltape at fas.harvard.edu>, "Anne Sieminski" <anne.sieminski at obs.ujfgrenoble.fr>, Paul Friberg <p.friberg at isti.com>, Kasper van Wijk <kasper at cgiss.boisestate.edu>, "Dylan Mikesell" <dmikesell at cgiss.boisestate.edu>, Federica Magnoni <federica.magnoni at ingv.it>, Ebru Bozdag <bozdag at princeton.edu>, Hejun Zhu <hejunzhu at princeton.edu>, Pieyre Le Loher <pieyre.le_loher at inria.fr>, Christina Morency <cmorency at Princeton.EDU>, Emanuele Casarotti <emanuele.casarotti at gmail.com>, Piero Basini <basini at bo.ingv.it>, "Emmanuel Chaljub" <Emmanuel.Chaljub at obs.ujfgrenoble.fr>, Qinya Liu <liuqy at physics.utoronto.ca>, Yang Luo <yangl at Princeton.EDU>, Daniel Peter <dpeter at Princeton.EDU>
+ fix the RK4 attenuation implementation problem from my 1999 approximation: will be done by Zhinan Xie
Hi all,
+ implement CPML: will be done by Zhinan Xie
Just to add that symplectic schemes are based on conservative systems
(they approximate the Hamiltonian) and we don't know yet how they behave
for dissipative media. But still of course something to be added and
looked at. Should also give credit to JeanPaul Ampuero who really first
suggested trying these schemes for elastodynamics with the SEM.
+ use PTScotch to be able to switch from 6 N^2 MPI threads to any other number, and feed that to SPECFEM3D instead of SPECFEM3D_GLOBE: will be done by Joseph Charles
Best regards,
Tarje
+=============================
On 23/07/10 02:21, Dimitri Komatitsch wrote:
>
> Hi Ying,
>
> What I know for sure if that back in 1999 when I developed the time
> scheme for attenuation I used a trick that made implementation much
> easier but that also makes the 4th order Runge Kutta (RK4) time scheme
> become second order only, i.e. RK2 instead of RK4. This means that for
> very long simulations (for instance multiorbit surface waves)
> attenuation can become inaccurate because of a lack of accuracy in the
> time scheme (in which case reducing Delta_t purposely solves the
> problem, but of course makes the simulation more expensive).
>
> A way of solving this problem could be to switch to better time
> schemes such as the symplectic time scheme that Tarje introduced in a
> GJI paper a few years ago. Tarje and I should probably implement that
> in the code one day... (therefore I cc him; we talked about this at
> the EGU meeting in May)
>
> Thank you,
> Dimitri.
+Less important for now:
+
 use RK4 or symplectic time scheme (will be done by Tarje NissenMeyer and/or JeanPaul Ampuero) (would be useful in the 2D version of the code as well):
Hi Jeroen, Perfect. I think talking to JeanPaul Ampuero would be useful
as well because in Utrecht last year he had told us that
he had implemented some nice 4thorder symplectic schemes
in his version of SEM2D. Dimitri.
Jeroen Tromp wrote:
> Hi Dimitri:
> This is one of the first things Tarje and I plan to work on after he
> arrives.
> Jeroen
> Dimitri Komatitsch wrote:
>> Hi Jeroen,
>> I think the last important thing that is missing
>> in SPECFEM3D (and SPECFEM2D) is a fourthorder time scheme.
>> I think both JeanPaul Ampuero and Tarje NissenMeyer
>> have worked on this, they will both be at Caltech
>> soon therefore maybe they could take care of adding it?
>> This would definitely increase the accuracy of very long
>> simulations (e.g. multiorbit surface waves).
>> Dimitri.

 From Dimitri: there is something in SPECFEM3D_GLOBE that we noticed a few years ago
regarding attenuation but never fixed: on page 813 of our 1999 paper
I used a trick suggested by Robertsson et al. (1994)
to use a nonstaggered RungeKutta (RK4) scheme for the attenuation
equations while using a staggered finitedifference (Newmark)
time scheme for all the other equations.
(combining staggered and nonstaggered formulations being difficult)
The trick works fine in most cases but there is a hidden problem for very long
simulations, and in particular for multiorbit surface waves:
because of that trick, the RK4 is not really a real RungeKutta scheme
(because grad(displacement) is used as a source but not known half way between
time steps, therefore I replaced it with an average between
t and t + Delta_t, which implies that the approximation is not fourthorder
accurate any more because of that smoothing/interpolation). Therefore
dispersion due to attenuation is not very accurately computed in the case of
very long runs; that matters mostly for surface waves.
We should fix that one day. The only thing to do would be to design a better time
integration scheme for the attenuation equation, without that trick; the rest
(Newmark etc) is fine and does not need to change.
Daniel Peter did not fix that in 2011; he only fixed the old implementation in SPECFEM3D which used a
fixed period absorption band to an adaptive one, which uses the resolved
periods of the mesh (as is done in the global version). The time stepping of the memory variables did not change.
Thus this IMPORTANT KNOWN BUG remains: in compute_forces_viscoelastic.f90 the calculation of
attenuation (viscoelasticity) is slightly incorrect because the gradients
are computed twice but *at the same time step* instead of at different
(staggered) time steps (t_{n1} and t_n or something like that).
That's easy to fix but I have no time for now. Let us fix that in the future.
It will only make a very small difference in the final seismograms (except for very long runs,
e.g., for multiorbit surface waves) therefore the current code with the bug can be used without any major problem.

Reply from Daniel Peter (Princeton University, USA) on February 23, 2010:
for higher orbit arrivals, we do might want to consider the symplectic time
schemes as well. Tarje made some tests and showed some "shocking" results
for longer distance paths. However, I am not sure if his scheme currently
considers purely elastic simulations, without incorporating attenuation.
It would be nice having both improved formulations included for those cases.
This seems to be a nice outlook for a future version update.

 use a better checkpointing/restarting system with CRC32 as in SPECFEM3D_GLOBE/tags/v4.1.0_beta_merged_mesher_solver_non_blocking_MPI/src
 use a potential of (rho * u) instead of u in the fluid, in case of fluidfluid discontinuities
 merge SPECFEM3D_GLOBE into SPECFEM3D, i.e. create a single package and see the mesh of the Earth as a general mesh that we would decompose in parallel using PTScotch. Problem: would this scale for huge global Earth simulations, e.g. at a seismic period of 1 second on 200,000 processor cores?
 make the code compatible with helioseismology / general Cowling formulation (could be done with Tarje NissenMeyer)
+ make the code compatible with helioseismology / general Cowling formulation
 could be done by Vala:
we could use heuristic rules to make source and receiver detection much
@@ 140,26 +56,6 @@


Things to do later:


 afficher le modele PREM,les courbes de nombre de points par longueur d'onde et de stabilite et les worst elements (stability et nb pts per lambda) en OpenDX en serial sur le maitre dans le nouveau mesher, en copiant/collant ce qu'il faut depuis display_prem_sampling_doubling.f90 sur le master seulement (rank == 0)






Things to add to the manual:


 ajouter les courbes de dispersion et stabilite au manuel pour tous les NEX classiques (160, 256, 380 etc) dans un Appendix special ladessus; dire dans cet appendix que dans le outer core ce n'est bien sur pas Vs qui est represente mais Vp / 1.25 (on a pris 25 % de marge car on travaille en velocity potential dans le outer core, ce qui demande un echantillonnage un peu meilleur)






Done:

More information about the CIGCOMMITS
mailing list