[cig-commits] r21737 - seismo/3D/SPECFEM3D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Apr 5 17:19:04 PDT 2013
Author: dkomati1
Date: 2013-04-05 17:19:04 -0700 (Fri, 05 Apr 2013)
New Revision: 21737
Modified:
seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
Log:
added missing GPU features to the to-do list
Modified: seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2013-04-05 22:22:33 UTC (rev 21736)
+++ seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2013-04-06 00:19:04 UTC (rev 21737)
@@ -45,7 +45,7 @@
the "elmnts_load" array up in decompose_mesh/part_decompose_mesh.f90 (around line 1432)
currently does not take into account C-PML elements weights. The matching array in the code is
-called "CPML_type" and is defined as follows: 1 = face, 2 = edge, 3 = corner.
+called "CPML_type" and is defined as follows: 1 = face, 2 = edge, 3 = corner.
that is similar to what Daniel does for acoustic elements in the current version of SPECFEM3D: he uses a weight of 1 for acoustic
elements and 3 (or something like that) for elastic elements, which use a 3D vector instead of a scalar.
@@ -79,7 +79,7 @@
That is the only file that would need to change (in addition to part_decompose_mesh.f90, which is called by decompose_mesh.F90).
Unfortunately that file is not very well written nor very well commented.
Thus unfortunately that set of routines is not so easy to understand and modify.
-A bunch of things are hardwired.
+A bunch of things are hardwired.
Thus, this is not so easy to do (Daniel and I checked about a year ago). It would maybe be easier to rewrite 'decompose_mesh_scotch'
from scratch as a new, MPI program instead of trying to add MPI support to the current and not so well written serial version.
@@ -122,12 +122,12 @@
several different types of sensitivity kernel files) into a single file. If so, some kernel processing scripts and tools will need to be adapted
accordingly.
-Feedback from Qinya: I can't agree more. We have had experience running big specfem jobs on our National Supercomputer Consortium
- (related to kernels for noise measurements, which generate even greater number of files). It crippled the entire file system,
+Feedback from Qinya: I can't agree more. We have had experience running big specfem jobs on our National Supercomputer Consortium
+ (related to kernels for noise measurements, which generate even greater number of files). It crippled the entire file system,
and the system administrators became really unfriendly to us after that ...
-On the other hand, I know a lot of visualization programs (paraview) actually read the individual binary files and combine them
- to form bigger visualization domains. So if we rewrite this, we need to write it in a form so that it is easy/faster for an
+On the other hand, I know a lot of visualization programs (paraview) actually read the individual binary files and combine them
+ to form bigger visualization domains. So if we rewrite this, we need to write it in a form so that it is easy/faster for an
external program to access bits of it (such as x,y,z, ibool, etc). Maybe direct-accessed C binary files?
@@ -188,7 +188,7 @@
- suggestion 30:
----------------
-Regarding the implicit time schemes that Zhinan has implemented in 2D, I agree that it would be great to put that in the official SVN version relatively soon to avoid losing the changes if we wait for too long. But I think we only need this in 2D for now, so let us not do it in 3D (at least for now in 2012). Let us just commit Zhinan's 2D version of the implicit routines to the official SVN code (making it off by default; the default should remain a purely explicit second-order Newmark scheme).
+Regarding the implicit time schemes that Zhinan has implemented in 2D, I agree that it would be great to put that in the official SVN version relatively soon to avoid losing the changes if we wait for too long. But I think we only need this in 2D for now, so let us not do it in 3D (at least for now in 2012). Let us just commit Zhinan's 2D version of the implicit routines to the official SVN code (making it off by default; the default should remain a purely explicit second-order Newmark scheme).
------------------------------------------------
@@ -221,16 +221,25 @@
Feedback from Qinya on the last point: I recently hosted a visit by Prof. Ling-yun Chiao (Chiao & Kuo 2001, Hung, Ping & Chiao 2011) from national Taiwan University, one of the first few people who started to apply wavelet to seismic inverse problem. I think we have some idea of how wavelet can potentially help with adjoint tomography, and I would love to share our thoughts with you all as well.
------------------------------------------------
-GPU:
+GPU version:
------------------------------------------------
- suggestion 34:
----------------
-DK DK new suggestion by Dimitri, Jan 2013: add GPU support for C-PML once the C-PML code is finished;
-C-PML is currently only implemented in regular Fortran + MPI
+Here are recent features of the code that have been implemented by different developers in the
+regular Fortran + MPI version but that have not been ported to the GPU version yet:
+- CPML
+- modified fluid/solid formulation of Yang Luo to allow for adjoint sources located in fluid regions
+ (for instance for offshore simulations in the oil industry)
+
+- QKappa attenuation support in addition to Qmu
+
+- new (fixed / improved / more accurate) attenuation time scheme of Zhinan Xie
+
+
- suggestion 22: (not really critical, probably not worth doing, at least for now)
----------------
@@ -268,7 +277,7 @@
this suggestion comes from notes that I took while talking to Paul Cristini a few months ago therefore I am not sure I remember
very well what Paul wanted to change about that. Paul, could you remind us more precisely the change that you had in mind?
-Answer from Paul:
+Answer from Paul:
My concern was to be able to have in one domain several media which can be elastic or visoelastic. Currently, there is one flag for the whole domain which makes all media elastic or visoelastic.
@@ -284,13 +293,13 @@
- suppress the TURN_ATTENUATION_ON flag
-- for each material in the list of materials, add a flag that says if that material is elastic or viscoelastic;
+- for each material in the list of materials, add a flag that says if that material is elastic or viscoelastic;
or detect if Q > 9000 for instance and if so, consider that that material is elastic and turn off attenuation for it
(the same way we detect acoustic media by checking if the S wave velocity Vs is smaller than 0.0001)
Otherwise, in the current implementation, we waste memory and CPU time by doing viscoelastic simulations everywhere, even in elastic regions, when TURN_ATTENUATION_ON is on.
-More precise analysis from Dimitri:
+More precise analysis from Dimitri:
this means in particular changing this in DATA/Par_file:
ATTENUATION_VISCOELASTIC_SOLID = .false. # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
@@ -317,25 +326,25 @@
------------------------------------------------
+ fault_solver_common:
- - make ordered version of subroutine assemble_MPI_vector_ext_mesh, and use it in subroutine initialize_fault
+ - make ordered version of subroutine assemble_MPI_vector_ext_mesh, and use it in subroutine initialize_fault
+ fault_solver_dynamic:
- - many hard-coded ad hoc features need to be set instead as user options
+ - many hard-coded ad hoc features need to be set instead as user options
+ rate-and-state friction:
- - make it a user-friendly option (currently a flag in fault_solver)
+ - make it a user-friendly option (currently a flag in fault_solver)
+ cubit interface:
- - consolidate python scripts (eliminate redundancy)
+ - consolidate python scripts (eliminate redundancy)
-+ verify consistency of fault edge nodes in kinematic solver:
++ verify consistency of fault edge nodes in kinematic solver:
currently, these nodes are not split but they are included in the fault database
+ simplify mesh generation with faults:
- - eliminate fault edge nodes?
+ - eliminate fault edge nodes?
Currently fault edge nodes must be defined as closed and non-split in CUBIT,
SESAME asigns a single ibulk, but it adds it to the fault database.
- - develop a recipe for mesh coarsening in CUBIT based on half "cubed sphere" mesh
+ - develop a recipe for mesh coarsening in CUBIT based on half "cubed sphere" mesh
- test tet-to-hex strategy
+ add snapshot outputs to kinematic solver (e.g. to export fault stresses for adjoint source inversion)
@@ -473,7 +482,7 @@
- make sure CPML can be used with external meshes (coming from Gmsh or CUBIT) as well, using the numbering convention
that you defined with Paul two months ago; make sure that external meshes with fluid regions and solid regions also work fine with CPML
-- clarify why the corners need additional arrays and additional code;
+- clarify why the corners need additional arrays and additional code;
I still do not understand that; I thought that corners were not different from the rest in PML and that memory variables were simply summed there,
but I see that in your code you have specific arrays for the corners, called "rmemory_dux_dx_corner()" and so on; they make the code
(the implementation) longer and more complex, thus let us see if they are absolutely needed or if we could (maybe?) get rid of them;
@@ -500,12 +509,12 @@
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
-See also if the current arrays called isPML() or PML() or something like that in SPECFEM3D need to be suppressed or on the contrary could
+See also if the current arrays called isPML() or PML() or something like that in SPECFEM3D need to be suppressed or on the contrary could
be reused. They seem to correspond to an earlier attempt a few years ago that was never finished.
-One thing that we will need to investigate in the specific case of GLOBE is how to make PML work for boundaries that are not aligned
-with x / y / z (i.e. for one chunk of SPECFEM3D_GLOBE). That should not be a problem but we will get more terms because of products with
-the three components of the normal vector. We already have the general (tensorial) formulation written in our PML paper of 2003,
+One thing that we will need to investigate in the specific case of GLOBE is how to make PML work for boundaries that are not aligned
+with x / y / z (i.e. for one chunk of SPECFEM3D_GLOBE). That should not be a problem but we will get more terms because of products with
+the three components of the normal vector. We already have the general (tensorial) formulation written in our PML paper of 2003,
but I never implemented all the terms.
The standard logical bit handling functions of Fortran are described for instance at http://www.ews.uiuc.edu/~mrgates2/docs/fortran.html
@@ -675,30 +684,30 @@
could you also rename "USE_RICKER_IPATI" to "USE_RICKER" or "USE_RICKER_TIME_FUNCTION" everywhere in the source code and in the input file?
Because "Ipati" is the name of an industrial field and thus we should make the name of that parameter less specific.
-Let us also move that parameter from the constants.h file to the input file (Par_file), to avoid having to recompile the code
+Let us also move that parameter from the constants.h file to the input file (Par_file), to avoid having to recompile the code
(see my previous email about this)
================
-As discussed yesterday during our Skype call, let us make the
+As discussed yesterday during our Skype call, let us make the
modifications below in the input Par_file of SPECFEM3D.
-NGNOD should be added and set to 8 or 27, and then when reading the
-Par_file the code should check that it is equal to one or the other, and
+NGNOD should be added and set to 8 or 27, and then when reading the
+Par_file the code should check that it is equal to one or the other, and
otherwise exit with an error message.
The code should then automatically set NGNOD2D to 4 if NGNOD == 8 and to 9 if NGNOD == 27.
-For the new movie flags, I suggest changing MOVIE_SURFACE and
-EXTERNAL_MESH_MOVIE_SURFACE, which are confusing, to MOVIE_SURFACE only,
-and then define a second (integer) parameter, called for instance
-MOVIE_TYPE, which could be equal to 1 to show the top surface
-(topography + oceans) only, to 2 to show all the external faces of the
+For the new movie flags, I suggest changing MOVIE_SURFACE and
+EXTERNAL_MESH_MOVIE_SURFACE, which are confusing, to MOVIE_SURFACE only,
+and then define a second (integer) parameter, called for instance
+MOVIE_TYPE, which could be equal to 1 to show the top surface
+(topography + oceans) only, to 2 to show all the external faces of the
mesh (i.e. topography + vertical edges + bottom), and so on.
Let me add this to the todo list. Let us also decide who does it
-(since new flags will be needed in the Par_file for PML, I guess it will
+(since new flags will be needed in the Par_file for PML, I guess it will
be easy and better to do all the modifications at the same time).
-Let us also modify the code to allow for shakemaps and movies at the
+Let us also modify the code to allow for shakemaps and movies at the
same time, rather than having to run the code twice.
> add NGNOD = 27 or 8 to the input Par_file (move it from shared/constants.h.in to the Par_file, since now either value can be used depending on the mesh,
@@ -728,7 +737,7 @@
- suggestion 05:
----------------
-Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
+Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
Date: Tue, 11 Sep 2012 18:34:21 +0200
From: Dimitri Komatitsch
To: Christina Morency
@@ -748,13 +757,13 @@
On 09/11/2012 06:17 PM, Christina Morency wrote:
> Hi Dimitri,
->
+>
> Yes, I can confirm what Yang is saying, it is already in the package,
> with interpolation if the source is not located on a GLL point (as it is
> indeed done in the 2D btw).
->
+>
> Christina
->
+>
> On 09/11/2012 09:14 AM, Yang Luo wrote:
>> I think Daniel (or someone else) has already copied necessary lines from old package.
>>
@@ -776,46 +785,46 @@
Feedback from Qinya: I think this is a great idea. In the past, we always had multiple versions of the same specfem code, just to be able to run point force sources. We can add a parameter `SOURCE_TYPE' in Par file, which is = 0 for moment-tensor source ( reads CMTSOLUTION, and assumes Heavside/err as source time function), and = 1 for point force source (and read something like PNTSOLUTION, including force location, direction, time shift and half duration). Now there is no standard stf to use for point force, but Heavside function does not seem to be a good idea because it will produce a static offset on all seismograms, which is not desirable. On the other hand, I don't see why Richer wavelet (2nd order derivative of Gaussian) is necessary better than Gaussian itself. One way to give user the freedom of choice is to put the stf info also in the PNTSOLUTION.
-Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
+Subject: Re: [CIG-SEISMO] simulation of seismic wave propagation excited by a force
Date: Tue, 11 Sep 2012 18:01:38 +0200
From: Dimitri Komatitsch
To: Yingzi Ying
Dear Yingzi,
-You are right, in the current code there is no option in the Par_file to
-use a point source, but we added this to the todo list right before the
-summer and thus we should/will definitely do it. It is easy to do: just
-add the Ricker wavelet to the accel() vector (which is in fact a force
-vector before we divide it by the mass matrix), to the component you
+You are right, in the current code there is no option in the Par_file to
+use a point source, but we added this to the todo list right before the
+summer and thus we should/will definitely do it. It is easy to do: just
+add the Ricker wavelet to the accel() vector (which is in fact a force
+vector before we divide it by the mass matrix), to the component you
want (vertical or horizontal) and you are all set.
-Joseph, could you maybe do it? below are a few lines from one of my test
-codes, in which I do just that; let us cut and paste these lines in
-SPECFEM3D (it is just a matter of selecting the right grid point; which
-we already do for a CMT source, thus we could add a force instead of a
-CMT at that location; the only difference being that the CMT source is
-added to all the points of a given spectral element, while a force
+Joseph, could you maybe do it? below are a few lines from one of my test
+codes, in which I do just that; let us cut and paste these lines in
+SPECFEM3D (it is just a matter of selecting the right grid point; which
+we already do for a CMT source, thus we could add a force instead of a
+CMT at that location; the only difference being that the CMT source is
+added to all the points of a given spectral element, while a force
vector should be added to a single point).
-Let me talk to Joseph tomorrow to see if we can release a patch of the
+Let me talk to Joseph tomorrow to see if we can release a patch of the
source code soon.
-I also cc Yang to see if he has already done that in a local version at
-Princeton maybe (for instance to model a force for some of his oil
+I also cc Yang to see if he has already done that in a local version at
+Princeton maybe (for instance to model a force for some of his oil
industry simulations in foothill regions?).
-Another difficulty is if the location of the force does not correspond
+Another difficulty is if the location of the force does not correspond
to a grid point; then there are two options:
-- select the closest grid point; i.e. slightly change the location of
+- select the closest grid point; i.e. slightly change the location of
the source (can be done automatically by the code)
-- use Christina Morency's way of adding a force source between grid
-points; I think I remember she did that in the 2D code, but I am not
+- use Christina Morency's way of adding a force source between grid
+points; I think I remember she did that in the 2D code, but I am not
sure how; thus I cc her.
-(that second option would be better, if Christina can confirm that it
+(that second option would be better, if Christina can confirm that it
works fine)
Best wishes,
@@ -863,15 +872,15 @@
Hi all,
Yes, let us definitely add it and make it an option in the Par_file.
-This will be very useful also in the case of non destructive testing of
+This will be very useful also in the case of non destructive testing of
materials (in which case the source is almost always a set of forces).
-What probably exists is a routine to add a force at a GLL grid point
-(which is easy to do). If we also want to have an option to model a
-source located somewhere between grid points, then I guess we will need
-to see how the discrete source term looks like (I think Daniel once told
-me that Christina had done it in the 2D code; probably summing the
-Lagrange interpolants at all the GLL points of the element? (as we do
+What probably exists is a routine to add a force at a GLL grid point
+(which is easy to do). If we also want to have an option to model a
+source located somewhere between grid points, then I guess we will need
+to see how the discrete source term looks like (I think Daniel once told
+me that Christina had done it in the 2D code; probably summing the
+Lagrange interpolants at all the GLL points of the element? (as we do
for receivers located between grid points?).
Thanks,
@@ -912,7 +921,7 @@
Hi Carl and Daniel,
Yes, I think that would be very useful;
-and we should also put these parameters in the input file (Par_file)
+and we should also put these parameters in the input file (Par_file)
rather than in constants.h, so that users do not need to recompile the code.
Thanks,
@@ -990,7 +999,7 @@
as Zhinan says, it is probably simple for Jo to make the modification by taking the code Zhinan added to modify
"rmass_x" and making the same modifications for "rmass_acoustic"
-(only when absorbing conditions are turned on of course).
+(only when absorbing conditions are turned on of course).
should be done in both SPECFEM3D and SPECFEM3D_GLOBE; already done in SPECFEM2D
@@ -1051,10 +1060,10 @@
this applies to SPECFEM3D_GLOBE; I am not sure about SPECFEM3D, please doublecheck
-Feedback from Qinya: I chose to output the exact nspec/nglob numbers for each slice to array_dims.txt in save_array_solver.txt,
- because the nspec/nglob number in values_from_mesher.h are only the maximum values for all slices. We need the exact numbers
- to read properly the topology/kernel binary files for combine_vol_data.f90, and generate the .mesh and .vtu files for paraview.
- If the new version of the code has exactly the same number of nspec/nglob number for all slices, it is no longer needed.
+Feedback from Qinya: I chose to output the exact nspec/nglob numbers for each slice to array_dims.txt in save_array_solver.txt,
+ because the nspec/nglob number in values_from_mesher.h are only the maximum values for all slices. We need the exact numbers
+ to read properly the topology/kernel binary files for combine_vol_data.f90, and generate the .mesh and .vtu files for paraview.
+ If the new version of the code has exactly the same number of nspec/nglob number for all slices, it is no longer needed.
Otherwise, I am sure there is a way around it. As long as we can get this info from other files, we can eliminate the array_dims.txt file.
------------------------------------------------
@@ -1134,8 +1143,8 @@
- suggestion 31:
----------------
-Elliott Sales de Andrade and Qinya will implement a routine to interpolate and output the sensitiviy kernels
-on a topologically-regular grid instead of the non-structured SEM grid (this can be useful for many processing
+Elliott Sales de Andrade and Qinya will implement a routine to interpolate and output the sensitiviy kernels
+on a topologically-regular grid instead of the non-structured SEM grid (this can be useful for many processing
packages that require a topologically-regular grid of points as input).
More information about the CIG-COMMITS
mailing list