[cig-commits] r20079 - seismo/3D/SPECFEM3D/trunk

[dkomati1 at geodynamics.org]

Author: dkomati1
Date: 2012-05-10 15:56:07 -0700 (Thu, 10 May 2012)
New Revision: 20079

Modified:
   seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
Log:
cleaned the old todo lists


Modified: seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt	2012-05-10 22:48:35 UTC (rev 20078)
+++ seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt	2012-05-10 22:56:07 UTC (rev 20079)
@@ -3,127 +3,51 @@
 ------------------------------------------------
 
 - suggestion 01:
-Zhinan commits his current clean 2D C-PML code to SVN; because changes 
-made outside of SVN are useless for regular users of the application, 
+Zhinan commits his current clean 2D C-PML code to SVN; because changes
+made outside of SVN are useless for regular users of the application,
 and are also quickly lost
 
 - suggestion 02:
-Max merges the GPU / CUDA improvements made by Peter Messmer into 
-SUNFLOWER; Daniel and Jo then merge SUNFLOWER into SPECFEM3D to get a 
-single version that supports both CPUs and GPUs (more precisely: either 
+Max merges the GPU / CUDA improvements made by Peter Messmer into
+SUNFLOWER; Daniel and Jo then merge SUNFLOWER into SPECFEM3D to get a
+single version that supports both CPUs and GPUs (more precisely: either
 CPUs or GPUs, depending on the case; no mixed / hybrid model);
 we then "svn delete" SUNFLOWER and keep a single code
 
-Older to-do list of SPECFEM3D (some items are probably still useful to implement):
-----------------------------------------------------------------------------------
 
-- 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 Cosserat or for some C-PML formulations this could matter
+- see (with Emanuele and also Paul Cristini) how to define and handle the PML absorbing elements in CUBIT; read this from input files in the solver
 
-- add C-PML
-
-- add an option to use a point source (a force) instead of a CMTSOLUTION, as suggested by Peng Chao from China.
-
-- automatic detection (coloring / flags) of the PML absorbing elements in CUBIT; read this from input files in the solver
-
-- Partitioning should be done using a parallel partitioner such as PT-SCOTCH.
-Developers and users in Switerland recently mentioned that using a serial domain decomposer such as SCOTCH was very slow for large meshes.
-Reading the (distributed or not) mesh, building the graph, distributing, and partitioning is not a problem, but redistributing the graph/mesh according to the partitioning obtained might well be. Should add support for PT-Scotch (i.e. decompose the mesh in parallel rather than serial; this will be required for very large meshes, too large to fit on a serial machine). Serial domain decomposition is already a limitation in the case of regional seismology, because serial meshes can be very large. One option would be to switch to ParMetis or PT-Scotch, i.e., the parallel version of these domain decomposition packages. I know that PT-Scotch works great and is pretty fast. Another option is ZOLTAN from Sandia National Labs, but I do not know if it is serial or parallel. I think adding support for PT-Scotch should not be too difficult if we give it a very simple initial "analytical" partition, e.g. if we have NSPEC spectral elements and NPROC domains, put the first 1..NPROC elements in domain 0, then NPROC+1...2*NPROC in domain 1 etc. That's your initial domain decomposition; then give that to PT-Scotch and ask it to optimize by moving elements between partitions.
-
-- there is something in SPECFEM3D 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 non-staggered Runge-Kutta (RK4) scheme for the attenuation
-equations while using a staggered finite-difference (Newmark)
-time scheme for all the other equations.
-(combining staggered and non-staggered formulations being difficult)
-The trick works fine in most cases but there is a hidden problem for very long
-simulations, and in particular for multi-orbit surface waves:
-because of that trick, the RK4 is not really a real Runge-Kutta 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 fourth-order
-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_{n-1} 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 therefore
-the current code with the bug can be used without any major problem.
-
-- add a checkpointing/restarting system with CRC-32 as in SPECFEM3D_GLOBE/tags/v4.1.0_beta_merged_mesher_solver_non_blocking_MPI/src
-
-- For binary files (for instance Carl Tape's m16 model) we should use NETCDF to store them and read them back; that is the only clean way of ensuring portability between different systems; and NETCDF is free and open source therefore there is no reason not to do it. 
-
-- we could add support for the TetGen mesh creation package, and cut each
-tetrahedron into four hexahedra using the barycenter. This could help design
-meshes for sedimentary basins, in particular near basin edges. Celine Blitz is probably going to work on that in Oct, Nov and Dec 2010.
-
-- we could add support for different polynomial degrees in different elements (p-adaptivity)
-
-- we could add support for Discontinuous Galerkin, and/or for SEM on tetrahedra and pyramids
-
-- Regarding memory size (getting an estimate of memory consumption in SPECFEM3D),
-somebody should just cut and paste my SPECFEM3D_GLOBE routine
-"SPECFEM3D_GLOBE/version41_beta/src/memory_eval.f90", which
-I call from "SPECFEM3D_GLOBE/version41_beta/src/create_header_file.f90".
-It would work for SPECFEM3D as well (with minor modifications).
-
-- add SOURCESOLUTION in addition to CMTSOLUTION, to choose between a CMT source and a force source
-
 - put a flag to choose between a Ricker source and a Heaviside source when a force source is used
 (we can cut and paste the Heaviside implementation from the case of a CMT source)
 
-- etudier comment lire un maillage CUBIT stocke au format natif de CUBIT (Exodus, base sur netCDF) depuis SPECFEM3D. On pourrait utiliser la commande "ncdump" dans CUBIT si necessaire d'apres Emanuele Casarotti. Une autre option serait d'utiliser le format de stockage ABAQUS dans CUBIT.
-
-
-by Nicolas Le Goff :
-
-Add-ONS :
-
- - Receivers distance to source along the surface computed for an asteroid can be done by computing the shortest path in a graph. No problem in serial, doing so in parallel might be tricky.
-
- - Distinguish between elements in contact with the envelope, and those that are in contact with a fracture : source, receivers, ... anything based on the envelope detection will have to take into account those fractures. Really tricky without prior information on the mesh.
-
- - Settle the issue of normal recording for the receivers. The normal is unique, but the plane can be described using two vectors at random (while conserving an orthonormal reference), so we have no reason to choose a pair of vectors compared to the others.i
-
-
-Older to-do list of SPECFEM3D_GLOBE (some items are probably still useful to implement):
-----------------------------------------------------------------------------------------
-
 - 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 --print-size --size-sort --reverse-sort --radix=d ./bin/xspecfem3D | cut -f 2- -d ' ' | less 
+nm --print-size --size-sort --reverse-sort --radix=d ./bin/xspecfem3D | cut -f 2- -d ' ' | less
 I think switching from (NGLLX,NGLLY,NGLLZ) to (1,1,1) i.e. using the lower left corner of each element in the case of 1D attenuation would be easy to implement and sufficient to fix the problem (since array sizes would decrease by a factor of 125).
 
-- 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 C-PML formulations this could matter
 
 - implement RK4 and LDDRK4-6 high-order time schemes, including for viscoelastic media: will be done by Zhinan Xie
 
-- fix the RK4 attenuation implementation problem from my 1999 approximation: will be done by Zhinan Xie
+Older to-do lists of SPECFEM2D, SPECFEM3D and SPECFEM3D_GLOBE (some items are probably still useful to implement):
+------------------------------------------------------------------------------------------------------------------
 
-- implement C-PML: will be done by Zhinan Xie
+(items listed in no particular order)
 
-- use PT-Scotch 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
+- add a checkpointing/restarting system with CRC-32 as in SPECFEM3D_GLOBE/tags/v4.1.0_beta_merged_mesher_solver_non_blocking_MPI/src , or using the fault tolerance library developed by Leonardo Bautista and Franck Cappello at INRIA
 
-=============================
+- Regarding memory size (getting an estimate of memory consumption in SPECFEM3D),
+somebody should just cut and paste my SPECFEM3D_GLOBE routine
+"SPECFEM3D_GLOBE/version41_beta/src/memory_eval.f90", which
+I call from "SPECFEM3D_GLOBE/version41_beta/src/create_header_file.f90".
+It would work for SPECFEM3D as well (with minor modifications).
 
-Less important for now:
-----------------------
+- For binary files (for instance Carl Tape's m16 model) we should use NETCDF to store them and read them back; that is the only clean way of ensuring portability between different systems; and NETCDF is free and open source therefore there is no reason not to do it.
 
-- use a better checkpointing/restarting system with CRC-32 as in SPECFEM3D_GLOBE/tags/v4.1.0_beta_merged_mesher_solver_non_blocking_MPI/src
+- in SPECFEM3D_GLOBE use a potential of (rho * u) instead of u in the fluid, in case of fluid-fluid discontinuities; probably already
+done in SPECFEM3D because of purely acoustic oil industry problems, but not done in SPECFEM3D_GLOBE because there is no first-order discontinuity inside
+the fluid outer core therefore it is not needed; thus this point is not needed except if we decide to use SPECFEM3D_GLOBE for other problems one day,
+for instance helioseismology
 
-- use a potential of (rho * u) instead of u in the fluid, in case of fluid-fluid 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 PT-Scotch. 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
@@ -146,48 +70,21 @@
 careful analysis we have found a way of doing better and getting rid of all
 the GATHERs and *all* the extra memory (this way we won't have to limit the number of sources to a maximum of 100000, as currently done in the solver).
 
-
----------------------------------------------------------------------------------
----------------------------------------------------------------------------------
----------------------------------------------------------------------------------
----------------------------------------------------------------------------------
-
-Done:
------
-
-- Jeroen Ritsema (from Univ of Michigan, USA) has an interesting suggestion
-(on February 22, 2010) for SPECFEM3D_GLOBE, an option that is currently missing
-but that would be easy to add, and potentially very useful:
-"If I may suggest one modification to specfem-3D: allow for 3D mantle
-structure (i.e. s20rts) and a PREM crust.
-This is helpful in isolating the mantle contribution to waveform anomalies.
-Right now, crust2.0 structure is automatically included if the mantle model is
-turned on." Solution implemented by Daniel Peter from Princeton Univ (USA) on February 23, 2010:
-"I put an option such that if one appends "_1Dcrust" to the 3D model
-name, e.g. "s20rts_1Dcrust" instead of "s20rts" in the Par_file, it
-will just take the crust from the corresponding 1D reference model
-(which is PREM for his s20rts model). Daniel. "
-
-- Purposely not done because we decided not to merge the mesher with the solver, after trying in version 4.1_beta and noticing that the merged code was difficult to write and to maintain: supprimer sections qui decrivent write_AVS_mesh_quality, check_buffers*.f90, check_mesh_quality*.f90 etc du manuel une fois que le mesher et le solver auront ete fusionnes
-
-
-
-Older to-do list of SPECFEM2D (some items are probably still useful to implement):
-----------------------------------------------------------------------------------
-
-- partitioning using weights for load-balancing
-
-- checking for points with different normals for absorbing conditions, when the absorbing edges are not in the same elements (similar to what is done for the corners).
-
-- SOMETHING THAT COULD BE MADE MORE GENERAL:
+- something that could be made more general in the 2D code:
 at line 770 of specfem2D.F90:
 !! DK DK if needed in the future, here the quality factor could be different for each point
 i.e. they could be given at each (i,j,ispec) instead of at each (ispec) only
 in the current version. Very easy to do if needed, just that line to change.
 
+- we could add support for the TetGen mesh creation package, and cut each
+tetrahedron into four hexahedra using the barycenter. This could help design
+meshes for sedimentary basins, in particular near basin edges. Celine Blitz is probably going to work on that in Oct, Nov and Dec 2010.
 
+- we could add support for different polynomial degrees in different elements (p-adaptivity)
 
+- we could add support for Discontinuous Galerkin, and/or for SEM on tetrahedra and pyramids
 
+
 ---------------------------------------------------------------------------------
 ---------------------------------------------------------------------------------
 ---------------------------------------------------------------------------------