[cig-commits] r22658 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: . EXAMPLES/small_benchmark_run_to_test_more_complex_Earth EXAMPLES/small_benchmark_run_to_test_very_simple_Earth doc/USER_MANUAL doc/USER_MANUAL/figures setup src/compute_optimized_dumping_undo_att src/shared src/specfem3D utils
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Jul 24 07:02:17 PDT 2013
Author: dkomati1
Date: 2013-07-24 07:02:17 -0700 (Wed, 24 Jul 2013)
New Revision: 22658
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_CSC_China.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Caltech.png
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_ETH.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Fairbanks.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Harvard.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_IPGP.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_IUF.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Intel_Exascale_Labs.png
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Maison_Simulation.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_NVIDIA.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_UPPA.png
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Univ_Toulouse.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_University_of_Toronto.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_aix_marseille_universite.png
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_cnrs.png
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_inria.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_princeton.pdf
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/expand_compute_strain_product.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/test_shape_for_vectorization_if_not_first_index.f90
Removed:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/Par_file
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/Par_file
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/specfem_3d_globe-cover_small.jpg
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/flags.guess
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/compute_optimized_dumping_undo_att/compute_optimized_dumping_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_classical.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_classical.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part3_kernel_computation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk
Log:
performed obvious merges from the trunk to the GPU branch
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/Par_file 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_more_complex_Earth/Par_file 2013-07-24 14:02:17 UTC (rev 22658)
@@ -57,10 +57,36 @@
# to undo attenuation for sensitivity kernel calculations or forward runs with SAVE_FORWARD
# use one (and only one) of the two flags below. UNDO_ATTENUATION is much better (it is exact)
# but requires a significant amount of disk space for temporary storage.
-PARTIAL_PHYS_DISPERSION_ONLY = .true.
+ATTENUATION_1D_WITH_3D_STORAGE = .true.
+PARTIAL_PHYS_DISPERSION_ONLY = .false.
UNDO_ATTENUATION = .false.
-NT_DUMP_ATTENUATION = 100 # how often we dump restart files to undo attenuation, only needed when using UNDO_ATTENUATION
+NT_DUMP_ATTENUATION = 77 # how often we dump restart files to undo attenuation, only needed when using UNDO_ATTENUATION
+# three mass matrices instead of one are needed to handle rotation very accurately;
+# otherwise rotation is handled slightly less accurately (but still reasonably well);
+# set to .true. if you are interested in precise effects related to rotation;
+# set to .false. if you are solving very large inverse problems at high frequency and also undoing attenuation exactly
+# using the UNDO_ATTENUATION flag above, in which case saving as much memory as possible can be a good idea.
+# You can also safely set it to .false. if you are not in a period range in which rotation matters, e.g. if you are targetting very short-period body waves.
+# if in doubt, set to .true.
+# Set it to .true. if you have ABSORBING_CONDITIONS above, because in that case the code will use the three mass matrices anyway
+# and thus there is no additional cost.
+# this flag is of course unused if ROTATION above is set to .false.
+EXACT_MASS_MATRIX_FOR_ROTATION = .false.
+
+# this for LDDRK high-order time scheme instead of Newmark
+USE_LDDRK = .false.
+
+# the maximum CFL of LDDRK is significantly higher than that of the Newmark scheme,
+# in a ratio that is theoretically 1.327 / 0.697 = 1.15 / 0.604 = 1.903 for a solid with Poisson's ratio = 0.25
+# and for a fluid (see the manual of the 2D code, SPECFEM2D, Tables 4.1 and 4.2, and that ratio does not
+# depend on whether we are in 2D or in 3D). However in practice a ratio of about 1.5 to 1.7 is often safer
+# (for instance for models with a large range of Poisson's ratio values).
+# Since the code computes the time step using the Newmark scheme, for LDDRK we will simply
+# multiply that time step by this ratio when LDDRK is on and when flag INCREASE_CFL_FOR_LDDRK is true.
+INCREASE_CFL_FOR_LDDRK = .true.
+RATIO_BY_WHICH_TO_INCREASE_IT = 1.5d0
+
# save AVS or OpenDX movies
#MOVIE_COARSE saves movie only at corners of elements (SURFACE OR VOLUME)
#MOVIE_COARSE does not work with create_movie_AVS_DX
@@ -125,6 +151,47 @@
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
-# output kernels on a regular grid instead of on the mesh points
+#-----------------------------------------------------------
+#
+# adjoint kernel outputs
+#
+#-----------------------------------------------------------
+
+# this parameter must be set to .true. to compute anisotropic kernels
+# in crust and mantle (related to the 21 Cij in geographical coordinates)
+# default is .false. to compute isotropic kernels (related to alpha and beta)
+ANISOTROPIC_KL = .false.
+
+# output only transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+# rather than fully anisotropic kernels when ANISOTROPIC_KL above is set to .true.
+# means to save radial anisotropic kernels, i.e., sensitivity kernels for beta_v, beta_h, etc.
+SAVE_TRANSVERSE_KL_ONLY = .false.
+
+# output approximate Hessian in crust mantle region.
+# means to save the preconditioning for gradients, they are cross correlations between forward and adjoint accelerations.
+APPROXIMATE_HESS_KL = .false.
+
+# forces transverse isotropy for all mantle elements
+# (default is to use transverse isotropy only between MOHO and 220)
+# means we allow radial anisotropy between the bottom of the crust to the bottom of the transition zone, i.e., 660~km depth.
+USE_FULL_TISO_MANTLE = .false.
+
+# output kernel mask to zero out source region
+# to remove large values near the sources in the sensitivity kernels
+SAVE_SOURCE_MASK = .false.
+
+# output kernels on a regular grid instead of on the GLL mesh points (a bit expensive)
SAVE_REGULAR_KL = .false.
+#-----------------------------------------------------------
+
+# set to true to use GPUs
+GPU_MODE = .false.
+
+# set to true to use the ADIOS library for I/Os
+ADIOS_ENABLED = .false.
+ADIOS_FOR_FORWARD_ARRAYS = .true.
+ADIOS_FOR_MPI_ARRAYS = .true.
+ADIOS_FOR_ARRAYS_SOLVER = .true.
+ADIOS_FOR_AVS_DX = .true.
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/Par_file 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/small_benchmark_run_to_test_very_simple_Earth/Par_file 2013-07-24 14:02:17 UTC (rev 22658)
@@ -57,10 +57,36 @@
# to undo attenuation for sensitivity kernel calculations or forward runs with SAVE_FORWARD
# use one (and only one) of the two flags below. UNDO_ATTENUATION is much better (it is exact)
# but requires a significant amount of disk space for temporary storage.
+ATTENUATION_1D_WITH_3D_STORAGE = .true.
PARTIAL_PHYS_DISPERSION_ONLY = .false.
UNDO_ATTENUATION = .false.
-NT_DUMP_ATTENUATION = 100 # how often we dump restart files to undo attenuation, only needed when using UNDO_ATTENUATION
+NT_DUMP_ATTENUATION = 77 # how often we dump restart files to undo attenuation, only needed when using UNDO_ATTENUATION
+# three mass matrices instead of one are needed to handle rotation very accurately;
+# otherwise rotation is handled slightly less accurately (but still reasonably well);
+# set to .true. if you are interested in precise effects related to rotation;
+# set to .false. if you are solving very large inverse problems at high frequency and also undoing attenuation exactly
+# using the UNDO_ATTENUATION flag above, in which case saving as much memory as possible can be a good idea.
+# You can also safely set it to .false. if you are not in a period range in which rotation matters, e.g. if you are targetting very short-period body waves.
+# if in doubt, set to .true.
+# Set it to .true. if you have ABSORBING_CONDITIONS above, because in that case the code will use the three mass matrices anyway
+# and thus there is no additional cost.
+# this flag is of course unused if ROTATION above is set to .false.
+EXACT_MASS_MATRIX_FOR_ROTATION = .false.
+
+# this for LDDRK high-order time scheme instead of Newmark
+USE_LDDRK = .false.
+
+# the maximum CFL of LDDRK is significantly higher than that of the Newmark scheme,
+# in a ratio that is theoretically 1.327 / 0.697 = 1.15 / 0.604 = 1.903 for a solid with Poisson's ratio = 0.25
+# and for a fluid (see the manual of the 2D code, SPECFEM2D, Tables 4.1 and 4.2, and that ratio does not
+# depend on whether we are in 2D or in 3D). However in practice a ratio of about 1.5 to 1.7 is often safer
+# (for instance for models with a large range of Poisson's ratio values).
+# Since the code computes the time step using the Newmark scheme, for LDDRK we will simply
+# multiply that time step by this ratio when LDDRK is on and when flag INCREASE_CFL_FOR_LDDRK is true.
+INCREASE_CFL_FOR_LDDRK = .true.
+RATIO_BY_WHICH_TO_INCREASE_IT = 1.5d0
+
# save AVS or OpenDX movies
#MOVIE_COARSE saves movie only at corners of elements (SURFACE OR VOLUME)
#MOVIE_COARSE does not work with create_movie_AVS_DX
@@ -125,6 +151,47 @@
# print source time function
PRINT_SOURCE_TIME_FUNCTION = .false.
-# output kernels on a regular grid instead of on the mesh points
+#-----------------------------------------------------------
+#
+# adjoint kernel outputs
+#
+#-----------------------------------------------------------
+
+# this parameter must be set to .true. to compute anisotropic kernels
+# in crust and mantle (related to the 21 Cij in geographical coordinates)
+# default is .false. to compute isotropic kernels (related to alpha and beta)
+ANISOTROPIC_KL = .false.
+
+# output only transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+# rather than fully anisotropic kernels when ANISOTROPIC_KL above is set to .true.
+# means to save radial anisotropic kernels, i.e., sensitivity kernels for beta_v, beta_h, etc.
+SAVE_TRANSVERSE_KL_ONLY = .false.
+
+# output approximate Hessian in crust mantle region.
+# means to save the preconditioning for gradients, they are cross correlations between forward and adjoint accelerations.
+APPROXIMATE_HESS_KL = .false.
+
+# forces transverse isotropy for all mantle elements
+# (default is to use transverse isotropy only between MOHO and 220)
+# means we allow radial anisotropy between the bottom of the crust to the bottom of the transition zone, i.e., 660~km depth.
+USE_FULL_TISO_MANTLE = .false.
+
+# output kernel mask to zero out source region
+# to remove large values near the sources in the sensitivity kernels
+SAVE_SOURCE_MASK = .false.
+
+# output kernels on a regular grid instead of on the GLL mesh points (a bit expensive)
SAVE_REGULAR_KL = .false.
+#-----------------------------------------------------------
+
+# set to true to use GPUs
+GPU_MODE = .false.
+
+# set to true to use the ADIOS library for I/Os
+ADIOS_ENABLED = .false.
+ADIOS_FOR_FORWARD_ARRAYS = .true.
+ADIOS_FOR_MPI_ARRAYS = .true.
+ADIOS_FOR_ARRAYS_SOLVER = .true.
+ADIOS_FOR_AVS_DX = .true.
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_CSC_China.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_CSC_China.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Caltech.png
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Caltech.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_ETH.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_ETH.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Fairbanks.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Fairbanks.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Harvard.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Harvard.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_IPGP.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_IPGP.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_IUF.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_IUF.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Intel_Exascale_Labs.png
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Intel_Exascale_Labs.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Maison_Simulation.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Maison_Simulation.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_NVIDIA.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_NVIDIA.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_UPPA.png
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_UPPA.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Univ_Toulouse.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_Univ_Toulouse.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_University_of_Toronto.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_University_of_Toronto.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_aix_marseille_universite.png
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_aix_marseille_universite.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_cnrs.png
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_cnrs.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_inria.jpg
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_inria.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_princeton.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/figures/logo_princeton.pdf
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/manual_SPECFEM3D_GLOBE.tex 2013-07-24 14:02:17 UTC (rev 22658)
@@ -94,7 +94,8 @@
\title{\thispagestyle{empty}\textbf{SPECFEM3D\_GLOBE}\\
\textbf{User Manual}}
%
-\author{$\copyright$ Princeton University (USA) and CNRS / INRIA / University of Pau (France)\\
+\author{$\copyright$ Princeton University (USA) and CNRS / University of Marseille (France),\\
+ETH Z\"urich (Switzerland),\\
Version 5.1}
\date{\noindent \today}
@@ -144,6 +145,38 @@
The cover graphic of the manual was created
by Santiago Lombeyda from Caltech's Center for Advanced Computing Research (CACR) \urlwithparentheses{http://www.cacr.caltech.edu}.\\
+\section*{Current and past main participants or main sponsors}
+
+\begin{figure}[htbp]
+\noindent \begin{centering}
+\includegraphics[width=0.125\textwidth]{figures/logo_cnrs}\vspace*{2truemm}
+\includegraphics[width=0.125\textwidth]{figures/logo_princeton}\vspace*{2truemm}
+\includegraphics[width=0.15\textwidth]{figures/logo_aix_marseille_universite}\vspace*{0.02truemm}
+\includegraphics[width=0.12\textwidth]{figures/logo_CSC_China}\vspace*{0.02truemm}
+\includegraphics[width=0.15\textwidth]{figures/logo_ETH}\vspace*{2truemm}
+\includegraphics[width=0.15\textwidth]{figures/logo_inria}\vspace*{2truemm}
+\includegraphics[width=0.14\textwidth]{figures/logo_UPPA}
+\par\end{centering}
+
+\noindent \begin{centering}
+\includegraphics[width=0.15\textwidth]{figures/logo_NVIDIA}\vspace*{2truemm}
+\includegraphics[width=0.125\textwidth]{figures/logo_IUF}\vspace*{2truemm}
+\includegraphics[width=0.135\textwidth]{figures/logo_Caltech}\vspace*{2truemm}
+\includegraphics[width=0.135\textwidth]{figures/logo_Harvard}\vspace*{2truemm}
+\includegraphics[width=0.135\textwidth]{figures/logo_IPGP}\vspace*{2truemm}
+\par\end{centering}
+
+\noindent \begin{centering}
+\includegraphics[width=0.112\textwidth]{figures/logo_University_of_Toronto}\vspace*{2truemm}
+%\includegraphics[width=0.119\textwidth]{figures/logo_INGV}\vspace*{2truemm} %% did not work on the 2D version
+\includegraphics[width=0.112\textwidth]{figures/logo_Univ_Toulouse}\vspace*{2truemm}
+%\includegraphics[width=0.112\textwidth]{figures/logo_TOTAL}\vspace*{2truemm} %% did not work on the 2D version
+\includegraphics[width=0.130\textwidth]{figures/logo_Fairbanks}\vspace*{2truemm}
+\includegraphics[width=0.140\textwidth]{figures/logo_Intel_Exascale_Labs}\vspace*{2truemm}
+\includegraphics[width=0.130\textwidth]{figures/logo_Maison_Simulation}
+\par\end{centering}
+\end{figure}
+
\newpage{}
\tableofcontents{}
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/doc/USER_MANUAL/specfem_3d_globe-cover_small.jpg
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/flags.guess
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/flags.guess 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/flags.guess 2013-07-24 14:02:17 UTC (rev 22658)
@@ -30,7 +30,7 @@
# parallel file systems like SFS 3.2 / Lustre 1.8. If omitted
# I/O throughput lingers at 2.5 MB/s, with it it can increase to ~44 MB/s
# However it does not make much of a difference on NFS mounted volumes or with SFS 3.1.1 / Lustre 1.6.7.1
- DEF_FFLAGS="-O3 -DFORCE_VECTORIZATION -check nobounds -xHost -ftz -assume buffered_io -assume byterecl -align sequence -vec-report0 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage" # -mcmodel=medium -shared-intel
+ DEF_FFLAGS="-O3 -DFORCE_VECTORIZATION -check nobounds -xHost -ftz -assume buffered_io -assume byterecl -align sequence -vec-report0 -std03 -diag-disable 6477 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage" # -mcmodel=medium -shared-intel
# useful for debugging...
# for debugging: change -O3 -DFORCE_VECTORIZATION -check nobounds to -check all -debug -g -O0 -fp-stack-check -traceback -ftrapuv
#
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2013-07-24 14:02:17 UTC (rev 22658)
@@ -5,8 +5,8 @@
!
! Main authors: Dimitri Komatitsch and Jeroen Tromp
! Princeton University, USA
-! and University of Pau / CNRS / INRIA, France
-! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! and CNRS / INRIA / University of Pau, France
+! (c) Princeton University and CNRS / INRIA / University of Pau
! April 2011
!
! This program is free software; you can redistribute it and/or modify
@@ -46,6 +46,12 @@
! set to SIZE_DOUBLE to run in double precision (increases memory size by 2)
integer, parameter :: CUSTOM_REAL = @CUSTOM_REAL@
+! this for non blocking assembly
+ integer, parameter :: ELEMENTS_NONBLOCKING_CM_IC = 1500
+ integer, parameter :: ELEMENTS_NONBLOCKING_OC = 3000
+
+!*********************************************************************************************************
+
! if files on a local path on each node are also seen as global with same path
! set to .true. typically on a machine with a common (shared) file system, e.g. LUSTRE, GPFS or NFS-mounted /home
! set to .false. typically on a cluster of nodes with local disks used to store the mesh (not very common these days)
@@ -234,7 +240,14 @@
! the mesher) and use them for the computation of boundary kernel (in the solver)
logical, parameter :: SAVE_BOUNDARY_MESH = .false.
-! mesh coloring:
+!*********************************************************************************************************
+! added these parameters for the GPU version of the solver with mesh coloring
+
+! sort outer elements first and then inner elements
+! in order to use non blocking MPI to overlap communications
+ logical, parameter :: SORT_MESH_INNER_OUTER = .true.
+ integer, parameter :: NGNOD_HEXAHEDRA = 8
+
! add mesh coloring for the GPU + MPI implementation
logical, parameter :: USE_MESH_COLORING_GPU = .false.
integer, parameter :: MAX_NUMBER_OF_COLORS = 1000
@@ -255,25 +268,6 @@
!!
!!-----------------------------------------------------------
-! this parameter must be set to .true. to compute anisotropic kernels
-! in crust and mantle (related to the 21 Cij in geographical coordinates)
-! default is .false. to compute isotropic kernels (related to alpha and beta)
- logical, parameter :: ANISOTROPIC_KL = .false.
-
-! output only transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
-! rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true.
- logical, parameter :: SAVE_TRANSVERSE_KL = .false.
-
-! output approximate hessian in crust mantle region
- logical, parameter :: APPROXIMATE_HESS_KL = .false.
-
-! output kernel mask to zero out source region
- logical, parameter :: SAVE_SOURCE_MASK = .false.
-
-! forces transverse isotropy for all mantle elements
-! (default is to use transverse isotropy only between MOHO and 220 )
- logical, parameter :: USE_FULL_TISO_MANTLE = .false.
-
! regular kernel parameters
character(len=*), parameter :: PATHNAME_KL_REG = 'DATA/kl_reg_grid.txt'
integer, parameter :: NM_KL_REG_LAYER = 100
@@ -283,12 +277,6 @@
real, parameter :: KL_REG_MIN_LAT = -90.0
real, parameter :: KL_REG_MAX_LAT = +90.0
-! forces attenuation arrays to be fully 3D and defined on all GLL points
-! (useful for more accurate attenuation profiles and adjoint inversions)
-! (if set to .false., 3D attenuation is only used for models with 3D attenuation distributions)
- logical, parameter :: USE_3D_ATTENUATION_ARRAYS = .false.
-
-
!!-----------------------------------------------------------
!!
!! time stamp information
@@ -361,6 +349,7 @@
! Deville routines optimized for NGLLX = NGLLY = NGLLZ = 5
integer, parameter :: m1 = NGLLX, m2 = NGLLX * NGLLY
+ integer, parameter :: NGLLCUBE = NGLLX * NGLLY * NGLLZ
! mid-points inside a GLL element
integer, parameter :: MIDX = (NGLLX+1)/2
@@ -372,6 +361,7 @@
! a few useful constants
double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
+ double precision, parameter :: ONE_HALF = HALF, ONE_FOURTH = 0.25d0, ONE_EIGHTH = 0.125d0, ONE_SIXTEENTH = 0.0625d0
real(kind=CUSTOM_REAL), parameter :: &
ONE_THIRD = 1._CUSTOM_REAL/3._CUSTOM_REAL, &
@@ -458,10 +448,10 @@
integer, parameter :: REFERENCE_MODEL_PREM = 1
integer, parameter :: REFERENCE_MODEL_IASP91 = 2
integer, parameter :: REFERENCE_MODEL_1066A = 3
- integer, parameter :: REFERENCE_MODEL_AK135 = 4
- integer, parameter :: REFERENCE_MODEL_1DREF = 5
+ integer, parameter :: REFERENCE_MODEL_AK135F_NO_MUD = 4
+ integer, parameter :: REFERENCE_MODEL_1DREF = 5
integer, parameter :: REFERENCE_MODEL_JP1D = 6
- integer, parameter :: REFERENCE_MODEL_SEA1D = 7
+ integer, parameter :: REFERENCE_MODEL_SEA1D = 7
! define flag for 3D Earth model
integer, parameter :: THREE_D_MODEL_S20RTS = 1
@@ -477,16 +467,6 @@
integer, parameter :: THREE_D_MODEL_S40RTS = 11
integer, parameter :: THREE_D_MODEL_GAPP2 = 12
-! define flag for regions of the global Earth for attenuation
- integer, parameter :: NUM_REGIONS_ATTENUATION = 5
-
- integer, parameter :: IREGION_ATTENUATION_INNER_CORE = 1
- integer, parameter :: IREGION_ATTENUATION_CMB_670 = 2
- integer, parameter :: IREGION_ATTENUATION_670_220 = 3
- integer, parameter :: IREGION_ATTENUATION_220_80 = 4
- integer, parameter :: IREGION_ATTENUATION_80_SURFACE = 5
- integer, parameter :: IREGION_ATTENUATION_UNDEFINED = 6
-
! number of standard linear solids for attenuation
integer, parameter :: N_SLS = 3
@@ -496,15 +476,11 @@
integer, parameter :: ATTENUATION_COMP_RESOLUTION = 1
integer, parameter :: ATTENUATION_COMP_MAXIMUM = 5000
-! for lookup table for attenuation every 100 m in radial direction of Earth model
- integer, parameter :: NRAD_ATTENUATION = 70000
- double precision, parameter :: TABLE_ATTENUATION = R_EARTH_KM * 10.0d0
-
! for determination of the attenuation period range
! if this is set to .true. then the hardcoded values will be used
! otherwise they are computed automatically from the Number of elements
! This *may* be a useful parameter for Benchmarking against older versions
- logical, parameter :: ATTENUATION_RANGE_PREDEFINED = .false.
+ logical, parameter :: ATTENUATION_RANGE_PREDEFINED = .false.
! flag for the four edges of each slice and for the bottom edge
integer, parameter :: XI_MIN = 1
@@ -525,12 +501,12 @@
! number of faces a given slice can share with other slices
! this is at most 2, except when there is only once slice per chunk
! in which case it is 4
- integer, parameter :: NUMFACES_SHARED = 4
+ integer, parameter :: NUMFACES_SHARED = 2 !!!!! DK DK removed support for one slice only, was 4
! number of corners a given slice can share with other slices
! this is at most 1, except when there is only once slice per chunk
! in which case it is 4
- integer, parameter :: NUMCORNERS_SHARED = 4
+ integer, parameter :: NUMCORNERS_SHARED = 1 !!!!!! DK DK removed support for one slice only, was 4
! number of slaves per corner
integer, parameter :: NUMSLAVES = 2
@@ -576,14 +552,38 @@
! for lookup table for gravity every 100 m in radial direction of Earth model
integer, parameter :: NRAD_GRAVITY = 70000
-!!!!!!!!!!!!!! parameters added for the thread-safe version of the code
-
!!-----------------------------------------------------------
!!
!! model constants
!!
!!-----------------------------------------------------------
+! number of layers in DATA/1066a/1066a.dat
+ integer, parameter :: NR_1066A = 160
+
+! number of layers in ak135-F
+ integer, parameter :: NR_AK135F_NO_MUD = 136
+
+! number of layers in DATA/s362ani/REF
+ integer, parameter :: NR_REF = 750
+
+! number of layers in DATA/Lebedev_sea99 1D model
+ integer, parameter :: NR_SEA1D = 163
+
+! three_d_mantle_model_constants
+ integer, parameter :: NK_20 = 20,NS_20 = 20,NS_40 = 40, ND = 1
+
+! heterogen_mantle_model_constants
+ integer, parameter :: N_R = 256,N_THETA = 256,N_PHI = 256
+
+! Japan 3D model (Zhao, 1994) constants
+ integer, parameter :: MPA=42,MRA=48,MHA=21,MPB=42,MRB=48,MHB=18
+ integer, parameter :: MKA=2101,MKB=2101
+
+! QRFSI12 constants
+ integer,parameter :: NKQ=8,MAXL_Q=12
+ integer,parameter :: NSQ=(MAXL_Q+1)**2,NDEPTHS_REFQ=913
+
! The meaningful range of Zhao et al. (1994) model is as follows:
! latitude : 32 - 45 N
! longitude: 130-145 E
@@ -595,6 +595,30 @@
double precision,parameter :: LON_MIN = 130.d0
double precision,parameter :: DEP_MAX = 500.d0
+! crustal_model_constants
+! crustal model parameters for crust2.0
+ integer, parameter :: NKEYS_CRUST = 359
+ integer, parameter :: NLAYERS_CRUST = 8
+ integer, parameter :: NCAP_CRUST = 180
+
+! General Crustmaps parameters
+ integer, parameter :: CRUSTMAP_RESOLUTION = 4 !means 1/4 degrees
+ integer, parameter :: NLAYERS_CRUSTMAP = 5
+
+! parameters for EPCRUST , from Molinari & Morelli model(2011)
+! latitude : 9.0N - 89.5N
+! longitude: 56.0W - 70.0E
+ character(len=*), parameter :: PATHNAME_EPCRUST = 'DATA/epcrust/EPcrust_0_5.txt'
+ integer, parameter :: EPCRUST_NLON = 253, EPCRUST_NLAT = 162, EPCRUST_NLAYER = 3
+ double precision, parameter :: EPCRUST_LON_MIN = -56.0d0
+ double precision, parameter :: EPCRUST_LON_MAX = 70.0d0
+ double precision, parameter :: EPCRUST_LAT_MIN = 9.0d0
+ double precision, parameter :: EPCRUST_LAT_MAX = 89.5d0
+ double precision, parameter :: EPCRUST_SAMPLE = 0.5d0
+ logical, parameter :: flag_smooth_epcrust = .true.
+ integer, parameter :: NTHETA_EP = 4, NPHI_EP = 20
+ double precision, parameter :: cap_degree_EP = 0.2d0
+
! parameters for GLL model (used for iterative model inversions)
character(len=*), parameter :: PATHNAME_GLL_modeldir = 'DATA/GLL/'
@@ -605,9 +629,6 @@
! integer, parameter :: GLL_REFERENCE_1D_MODEL = REFERENCE_MODEL_1DREF
! integer, parameter :: GLL_REFERENCE_MODEL = THREE_D_MODEL_S362ANI
-
-!!!!!!!!!!!!!! end of parameters added for the thread-safe version of the code
-
!!-----------------------------------------------------------
!!
!! crustal stretching
@@ -673,20 +694,9 @@
integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
-! this for LDDRK
- logical, parameter :: USE_LDDRK = .false. ! .true.
+! for LDDRK high-order time scheme
+ integer, parameter :: NSTAGE = 6
-! the maximum CFL of LDDRK is significantly higher than that of the Newmark scheme,
-! in a ratio that is theoretically 1.327 / 0.697 = 1.15 / 0.604 = 1.903 for a solid with Poisson's ratio = 0.25
-! and for a fluid (see the manual of the 2D code, SPECFEM2D, Tables 4.1 and 4.2, and that ratio does not
-! depend on whether we are in 2D or in 3D). However in practice a ratio of about 1.5 to 1.7 is often safer
-! (for instance for models with a large range of Poisson's ratio values).
-! Since the SPECFEM3D_GLOBE code computes the time step using the Newmark scheme, for LDDRK we will simply
-! multiply that time step by this ratio when LDDRK is on and when flag INCREASE_CFL_FOR_LDDRK is true.
- logical, parameter :: INCREASE_CFL_FOR_LDDRK = .true.
- double precision, parameter :: RATIO_BY_WHICH_TO_INCREASE_IT = 1.5d0
-
- integer, parameter :: NSTAGE = 6
real(kind=CUSTOM_REAL), parameter, dimension(NSTAGE) :: ALPHA_LDDRK = &
(/0.0_CUSTOM_REAL,-0.737101392796_CUSTOM_REAL, -1.634740794341_CUSTOM_REAL,&
-0.744739003780_CUSTOM_REAL,-1.469897351522_CUSTOM_REAL,-2.813971388035_CUSTOM_REAL/)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/compute_optimized_dumping_undo_att/compute_optimized_dumping_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/compute_optimized_dumping_undo_att/compute_optimized_dumping_undo_att.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/compute_optimized_dumping_undo_att/compute_optimized_dumping_undo_att.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -43,13 +43,15 @@
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NSOURCES,NTSTEP_BETWEEN_FRAMES, &
NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY, &
+ ATT1,ATT2,ATT3,ATT4,ATT5
double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
- MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,RMOHO_FICTITIOUS_IN_MESHER
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,RMOHO_FICTITIOUS_IN_MESHER, &
+ RATIO_BY_WHICH_TO_INCREASE_IT
logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE, &
@@ -59,7 +61,11 @@
ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD, &
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,HONOR_1D_SPHERICAL_MOHO,WRITE_SEISMOGRAMS_BY_MASTER,&
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,SAVE_REGULAR_KL
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,SAVE_REGULAR_KL, &
+ USE_LDDRK,INCREASE_CFL_FOR_LDDRK,ANISOTROPIC_KL,SAVE_TRANSVERSE_KL_ONLY,APPROXIMATE_HESS_KL, &
+ USE_FULL_TISO_MANTLE,SAVE_SOURCE_MASK,GPU_MODE,ADIOS_ENABLED,ADIOS_FOR_FORWARD_ARRAYS, &
+ ADIOS_FOR_MPI_ARRAYS,ADIOS_FOR_ARRAYS_SOLVER,ADIOS_FOR_AVS_DX, &
+ EXACT_MASS_MATRIX_FOR_ROTATION,ATTENUATION_1D_WITH_3D_STORAGE
character(len=150) LOCAL_PATH,MODEL
@@ -84,7 +90,7 @@
double precision :: static_memory_size
integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUATION, &
NSPEC_INNER_CORE_ATTENUATION, &
NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
@@ -141,7 +147,11 @@
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY,&
- SAVE_REGULAR_KL,PARTIAL_PHYS_DISPERSION_ONLY,UNDO_ATTENUATION,NT_DUMP_ATTENUATION)
+ SAVE_REGULAR_KL,PARTIAL_PHYS_DISPERSION_ONLY,UNDO_ATTENUATION,NT_DUMP_ATTENUATION, &
+ USE_LDDRK,INCREASE_CFL_FOR_LDDRK,ANISOTROPIC_KL,SAVE_TRANSVERSE_KL_ONLY,APPROXIMATE_HESS_KL, &
+ USE_FULL_TISO_MANTLE,SAVE_SOURCE_MASK,GPU_MODE,ADIOS_ENABLED,ADIOS_FOR_FORWARD_ARRAYS, &
+ ADIOS_FOR_MPI_ARRAYS,ADIOS_FOR_ARRAYS_SOLVER,ADIOS_FOR_AVS_DX,RATIO_BY_WHICH_TO_INCREASE_IT, &
+ ATT1,ATT2,ATT3,ATT4,ATT5,EXACT_MASS_MATRIX_FOR_ROTATION,ATTENUATION_1D_WITH_3D_STORAGE)
! optimal dumping interval calculation can only be done when SIMULATION_TYPE == 3 in the Par_file,
! thus set it to that value here in this serial code even if it has a different value in the Par_file
@@ -166,7 +176,7 @@
ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUATION, &
NSPEC_INNER_CORE_ATTENUATION, &
NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
@@ -176,7 +186,9 @@
NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,ATT1,ATT2,ATT3, &
+ APPROXIMATE_HESS_KL,ANISOTROPIC_KL,NOISE_TOMOGRAPHY, &
+ NCHUNKS,USE_LDDRK,EXACT_MASS_MATRIX_FOR_ROTATION,static_memory_size)
NGLOB1D_RADIAL_TEMP(:) = &
(/maxval(NGLOB1D_RADIAL_CORNER(1,:)),maxval(NGLOB1D_RADIAL_CORNER(2,:)),maxval(NGLOB1D_RADIAL_CORNER(3,:))/)
@@ -222,6 +234,8 @@
print *,'How much memory (in GB) is installed on your machine per CPU core?'
print *,' (or per GPU card or per INTEL MIC Phi board)?'
print *,' (beware, this value MUST be given per core, i.e. per MPI thread, i.e. per MPI rank, NOT per node)'
+ print *,' (this value is for instance 4 GB on Tiger at Princeton, 2 GB on the non-GPU part of Titan at ORNL i.e. when using'
+ print *,' CPUs only there, 2 GB also on the machine used by Christina Morency, and 1.5 GB on the GPU cluster in Marseille)'
read(*,*) gigabytes_avail_per_core
if(gigabytes_avail_per_core < 0.1d0) stop 'less than 100 MB per core does not seem realistic; exiting...'
@@ -230,11 +244,11 @@
print *
print *,'What percentage of this total do you allow us to use, keeping in mind that you'
print *,'need to leave some memory available for the GNU/Linux system to run?'
- print *,' (a typical value is 90%; 92% to 95% is probably OK too; 85% is very safe)'
+ print *,' (a typical value is 90%; 92% is probably OK too; 85% is very safe)'
read(*,*) percentage_to_use_per_core
if(percentage_to_use_per_core < 50.d0) stop 'less than 50% does not seem realistic; exiting...'
- if(percentage_to_use_per_core > 96.d0) stop 'more than 96% is risky; exiting...'
+ if(percentage_to_use_per_core > 95.d0) stop 'more than 95% is risky; exiting...'
what_we_can_use_in_GB = gigabytes_avail_per_core * percentage_to_use_per_core / 100.d0
@@ -267,7 +281,8 @@
size_to_store_at_each_time_step = size_to_store_at_each_time_step / 1.d9
print *
- print *,'each time step to store in memory to undo attenuation requires storing ',size_to_store_at_each_time_step,' GB per core'
+ print *,'each time step to store in memory to undo attenuation will require storing ', &
+ size_to_store_at_each_time_step,' GB per core'
print *
print *,'*******************************************************************************'
@@ -298,7 +313,7 @@
if (ATTENUATION) then
! R_memory_crust_mantle
disk_size_of_each_dumping = disk_size_of_each_dumping + 5.d0*dble(N_SLS)*dble(NGLLX)* &
- dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
+ dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUATION*dble(CUSTOM_REAL)
! R_memory_inner_core
disk_size_of_each_dumping = disk_size_of_each_dumping + 5.d0*dble(N_SLS)*dble(NGLLX)* &
@@ -315,18 +330,18 @@
print *,'we will need to save a total of ',number_of_dumpings_to_do,' dumpings (restart files) to disk'
print *
- print *,'each dumping on the disk to undo attenuation requires storing ',disk_size_of_each_dumping,' GB per core'
+ print *,'each dumping on the disk to undo attenuation will require storing ',disk_size_of_each_dumping,' GB per core'
print *
- print *,'each dumping on the disk requires storing ',disk_size_of_each_dumping*NPROCTOT, &
+ print *,'each dumping on the disk will require storing ',disk_size_of_each_dumping*NPROCTOT, &
' GB for all cores'
print *
- print *,'ALL dumpings on the disk require storing ',disk_size_of_each_dumping*number_of_dumpings_to_do,' GB per core'
+ print *,'ALL dumpings on the disk will require storing ',disk_size_of_each_dumping*number_of_dumpings_to_do,' GB per core'
print *
print *,'*******************************************************************************'
- print *,'ALL dumpings on the disk require storing ', &
+ print *,'ALL dumpings on the disk will require storing ', &
disk_size_of_each_dumping*number_of_dumpings_to_do*NPROCTOT,' GB for all cores'
print *,' i.e. ',disk_size_of_each_dumping*number_of_dumpings_to_do*NPROCTOT/1000.d0,' TB'
print *,'*******************************************************************************'
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk 2013-07-24 14:02:17 UTC (rev 22658)
@@ -33,7 +33,7 @@
shared_OBJECTS = \
$O/auto_ner.o \
- $O/broadcast_compute_parameters.o \
+ $O/broadcast_computed_parameters.o \
$O/calendar.o \
$O/count_number_of_sources.o \
$O/create_name_database.o \
@@ -79,8 +79,8 @@
$O/auto_ner.o: ${SETUP}/constants.h $S/auto_ner.f90
${FCCOMPILE_CHECK} -c -o $O/auto_ner.o ${FCFLAGS_f90} $S/auto_ner.f90
-$O/broadcast_compute_parameters.o: ${SETUP}/constants.h $S/broadcast_compute_parameters.f90
- ${MPIFCCOMPILE_CHECK} -c -o $O/broadcast_compute_parameters.o ${FCFLAGS_f90} $S/broadcast_compute_parameters.f90
+$O/broadcast_computed_parameters.o: ${SETUP}/constants.h $S/broadcast_computed_parameters.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/broadcast_computed_parameters.o ${FCFLAGS_f90} $S/broadcast_computed_parameters.f90
$O/calendar.o: $S/calendar.f90
${FCCOMPILE_CHECK} -c -o $O/calendar.o ${FCFLAGS_f90} $S/calendar.f90
@@ -149,8 +149,8 @@
$O/rthetaphi_xyz.o: ${SETUP}/constants.h $S/rthetaphi_xyz.f90
${FCCOMPILE_CHECK} -c -o $O/rthetaphi_xyz.o ${FCFLAGS_f90} $S/rthetaphi_xyz.f90
-$O/save_header_file.o: ${SETUP}/constants.h $S/save_header_file.f90
- ${FCCOMPILE_CHECK} -c -o $O/save_header_file.o ${FCFLAGS_f90} $S/save_header_file.f90
+$O/save_header_file.o: ${SETUP}/constants.h $S/save_header_file.F90
+ ${FCCOMPILE_CHECK} -c -o $O/save_header_file.o ${FCFLAGS_f90} $S/save_header_file.F90
$O/spline_routines.o: ${SETUP}/constants.h $S/spline_routines.f90
${FCCOMPILE_CHECK} -c -o $O/spline_routines.o ${FCFLAGS_f90} $S/spline_routines.f90
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90 (from rev 22656, seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -0,0 +1,544 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and CNRS / INRIA / University of Pau, France
+! (c) Princeton University and CNRS / INRIA / University of Pau
+! April 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! save header file OUTPUT_FILES/values_from_mesher.h
+
+ subroutine save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ ELLIPTICITY,GRAVITY,ROTATION, &
+ OCEANS,ATTENUATION,ATTENUATION_NEW,ATTENUATION_3D, &
+ ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
+ INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP,&
+ static_memory_size, &
+ NSPEC2D_TOP,NSPEC2D_BOTTOM, &
+ NSPEC2DMAX_YMIN_YMAX,NSPEC2DMAX_XMIN_XMAX, &
+ NPROC_XI,NPROC_ETA, &
+ NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
+ SIMULATION_TYPE,SAVE_FORWARD,MOVIE_VOLUME)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, dimension(MAX_NUM_REGIONS) :: NSPEC, nglob
+
+ integer NEX_XI,NEX_ETA,NPROC,NPROCTOT,NCHUNKS,NSOURCES,NSTEP
+
+ logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_NEW,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE
+
+ double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
+ CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH
+
+ ! static memory size needed by the solver
+ double precision :: static_memory_size
+
+ integer, dimension(MAX_NUM_REGIONS) :: &
+ NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_YMIN_YMAX,NSPEC2DMAX_XMIN_XMAX
+ integer :: NPROC_XI,NPROC_ETA
+
+ integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
+ NSPEC2D_MOHO, NSPEC2D_400, NSPEC2D_670, NSPEC2D_CMB, NSPEC2D_ICB
+
+ integer :: SIMULATION_TYPE
+ logical :: SAVE_FORWARD,MOVIE_VOLUME
+
+ ! local parameters
+ double precision :: subtract_central_cube_elems,subtract_central_cube_points
+ ! for regional code
+ double precision x,y,gamma,rgt,xi,eta
+ double precision x_top,y_top,z_top
+ double precision ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD
+ ! rotation matrix from Euler angles
+ integer i,j,ix,iy,icorner
+ double precision rotation_matrix(3,3)
+ double precision vector_ori(3),vector_rotated(3)
+ double precision r_corner,theta_corner,phi_corner,lat,long,colat_corner
+ integer :: att1,att2,att3,att4,att5
+ integer :: ier
+ character(len=150) HEADER_FILE
+
+ ! copy number of elements and points in an include file for the solver
+ call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
+ open(unit=IOUT,file=HEADER_FILE,status='unknown',iostat=ier)
+ if( ier /= 0 ) stop 'error opening OUTPUT_FILES/values_from_mesher.h'
+
+ write(IOUT,*)
+
+ write(IOUT,*) '!'
+ write(IOUT,*) '! this is the parameter file for static compilation of the solver'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! mesh statistics:'
+ write(IOUT,*) '! ---------------'
+ write(IOUT,*) '!'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! number of chunks = ',NCHUNKS
+ write(IOUT,*) '!'
+
+! the central cube is counted 6 times, therefore remove 5 times
+ if(INCLUDE_CENTRAL_CUBE) then
+ write(IOUT,*) '! these statistics include the central cube'
+ subtract_central_cube_elems = 5.d0 * dble((NEX_XI/8))**3
+ subtract_central_cube_points = 5.d0 * (dble(NEX_XI/8)*dble(NGLLX-1)+1.d0)**3
+ else
+ write(IOUT,*) '! these statistics do not include the central cube'
+ subtract_central_cube_elems = 0.d0
+ subtract_central_cube_points = 0.d0
+ endif
+
+ write(IOUT,*) '!'
+ write(IOUT,*) '! number of processors = ',NPROCTOT ! should be = NPROC
+ write(IOUT,*) '!'
+ write(IOUT,*) '! maximum number of points per region = ',nglob(IREGION_CRUST_MANTLE)
+ write(IOUT,*) '!'
+! use fused loops on NEC SX
+ write(IOUT,*) '! on NEC SX, make sure "loopcnt=" parameter'
+ write(IOUT,*) '! in Makefile is greater than max vector length = ',nglob(IREGION_CRUST_MANTLE)*NDIM
+ write(IOUT,*) '!'
+
+ write(IOUT,*) '! total elements per slice = ',sum(NSPEC)
+ write(IOUT,*) '! total points per slice = ',sum(nglob)
+ write(IOUT,*) '!'
+
+ write(IOUT,'(1x,a,i1,a)') '! total for full ',NCHUNKS,'-chunk mesh:'
+ write(IOUT,*) '! ---------------------------'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! exact total number of spectral elements in entire mesh = '
+ write(IOUT,*) '! ',dble(NCHUNKS)*dble(NPROC)*dble(sum(NSPEC)) - subtract_central_cube_elems
+ write(IOUT,*) '! approximate total number of points in entire mesh = '
+ write(IOUT,*) '! ',dble(NCHUNKS)*dble(NPROC)*dble(sum(nglob)) - subtract_central_cube_points
+! there are 3 DOFs in solid regions, but only 1 in fluid outer core
+ write(IOUT,*) '! approximate total number of degrees of freedom in entire mesh = '
+ write(IOUT,*) '! ',dble(NCHUNKS)*dble(NPROC)*(3.d0*(dble(sum(nglob))) &
+ - 2.d0*dble(nglob(IREGION_OUTER_CORE))) &
+ - 3.d0*subtract_central_cube_points
+ write(IOUT,*) '!'
+
+! display location of chunk if regional run
+ if(NCHUNKS /= 6) then
+
+ write(IOUT,*) '! position of the mesh chunk at the surface:'
+ write(IOUT,*) '! -----------------------------------------'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! angular size in first direction in degrees = ',sngl(ANGULAR_WIDTH_XI_IN_DEGREES)
+ write(IOUT,*) '! angular size in second direction in degrees = ',sngl(ANGULAR_WIDTH_ETA_IN_DEGREES)
+ write(IOUT,*) '!'
+ write(IOUT,*) '! longitude of center in degrees = ',sngl(CENTER_LONGITUDE_IN_DEGREES)
+ write(IOUT,*) '! latitude of center in degrees = ',sngl(CENTER_LATITUDE_IN_DEGREES)
+ write(IOUT,*) '!'
+ write(IOUT,*) '! angle of rotation of the first chunk = ',sngl(GAMMA_ROTATION_AZIMUTH)
+
+! convert width to radians
+ ANGULAR_WIDTH_XI_RAD = ANGULAR_WIDTH_XI_IN_DEGREES * DEGREES_TO_RADIANS
+ ANGULAR_WIDTH_ETA_RAD = ANGULAR_WIDTH_ETA_IN_DEGREES * DEGREES_TO_RADIANS
+
+! compute rotation matrix from Euler angles
+ call euler_angles(rotation_matrix,CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH)
+
+! loop on the four corners of the chunk to display their coordinates
+ icorner = 0
+ do iy = 0,1
+ do ix = 0,1
+
+ icorner = icorner + 1
+
+ xi= - ANGULAR_WIDTH_XI_RAD/2. + dble(ix)*ANGULAR_WIDTH_XI_RAD
+ eta= - ANGULAR_WIDTH_ETA_RAD/2. + dble(iy)*ANGULAR_WIDTH_ETA_RAD
+
+ x=dtan(xi)
+ y=dtan(eta)
+
+ gamma=ONE/dsqrt(ONE+x*x+y*y)
+ rgt=R_UNIT_SPHERE*gamma
+
+! define the mesh points at the top surface
+ x_top=-y*rgt
+ y_top=x*rgt
+ z_top=rgt
+
+! rotate top
+ vector_ori(1) = x_top
+ vector_ori(2) = y_top
+ vector_ori(3) = z_top
+ do i=1,3
+ vector_rotated(i)=0.0d0
+ do j=1,3
+ vector_rotated(i)=vector_rotated(i)+rotation_matrix(i,j)*vector_ori(j)
+ enddo
+ enddo
+ x_top = vector_rotated(1)
+ y_top = vector_rotated(2)
+ z_top = vector_rotated(3)
+
+! convert to latitude and longitude
+ call xyz_2_rthetaphi_dble(x_top,y_top,z_top,r_corner,theta_corner,phi_corner)
+ call reduce(theta_corner,phi_corner)
+
+! convert geocentric to geographic colatitude
+ colat_corner=PI_OVER_TWO-datan(1.006760466d0*dcos(theta_corner)/dmax1(TINYVAL,dsin(theta_corner)))
+ if(phi_corner>PI) phi_corner=phi_corner-TWO_PI
+
+! compute real position of the source
+ lat = (PI_OVER_TWO-colat_corner)*RADIANS_TO_DEGREES
+ long = phi_corner*RADIANS_TO_DEGREES
+
+ write(IOUT,*) '!'
+ write(IOUT,*) '! corner ',icorner
+ write(IOUT,*) '! longitude in degrees = ',long
+ write(IOUT,*) '! latitude in degrees = ',lat
+
+ enddo
+ enddo
+
+ write(IOUT,*) '!'
+
+ endif ! regional chunk
+
+ write(IOUT,*) '! resolution of the mesh at the surface:'
+ write(IOUT,*) '! -------------------------------------'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! spectral elements along a great circle = ',4*NEX_XI
+ write(IOUT,*) '! GLL points along a great circle = ',4*NEX_XI*(NGLLX-1)
+ write(IOUT,*) '! average distance between points in degrees = ',360./real(4*NEX_XI*(NGLLX-1))
+ write(IOUT,*) '! average distance between points in km = ',real(TWO_PI*R_EARTH/1000.d0)/real(4*NEX_XI*(NGLLX-1))
+ write(IOUT,*) '! average size of a spectral element in km = ',real(TWO_PI*R_EARTH/1000.d0)/real(4*NEX_XI)
+ write(IOUT,*) '!'
+ write(IOUT,*) '! number of time steps = ',NSTEP
+ write(IOUT,*) '!'
+ write(IOUT,*) '! number of seismic sources = ',NSOURCES
+ write(IOUT,*) '!'
+ write(IOUT,*)
+
+ write(IOUT,*) '! approximate static memory needed by the solver:'
+ write(IOUT,*) '! ----------------------------------------------'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! (lower bound, usually the real amount used is 5% to 10% higher)'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! (you can get a more precise estimate of the size used per MPI process'
+ write(IOUT,*) '! by typing "size -d bin/xspecfem3D"'
+ write(IOUT,*) '! after compiling the code with the DATA/Par_file you plan to use)'
+ write(IOUT,*) '!'
+ write(IOUT,*) '! size of static arrays per slice = ',static_memory_size/1.d6,' MB'
+ write(IOUT,*) '! = ',static_memory_size/1048576.d0,' MiB'
+ write(IOUT,*) '! = ',static_memory_size/1.d9,' GB'
+ write(IOUT,*) '! = ',static_memory_size/1073741824.d0,' GiB'
+ write(IOUT,*) '!'
+
+ ! note: using less memory becomes an issue only if the strong scaling of the code is poor.
+ ! Some users will run simulations with an executable using far less than 80% RAM per core
+ ! if they prefer having a faster computational time (and use a higher number of cores).
+
+ write(IOUT,*) '! (should be below to 80% or 90% of the memory installed per core)'
+ write(IOUT,*) '! (if significantly more, the job will not run by lack of memory )'
+ write(IOUT,*) '! (note that if significantly less, you waste a significant amount'
+ write(IOUT,*) '! of memory per processor core)'
+ write(IOUT,*) '! (but that can be perfectly acceptable if you can afford it and'
+ write(IOUT,*) '! want faster results by using more cores)'
+ write(IOUT,*) '!'
+ if(static_memory_size*dble(NPROCTOT)/1.d6 < 10000.d0) then
+ write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1.d6,' MB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1048576.d0,' MiB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1.d9,' GB'
+ else
+ write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1.d9,' GB'
+ endif
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GiB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1.d12,' TB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TiB'
+ write(IOUT,*) '!'
+
+ write(IOUT,*)
+ write(IOUT,*) 'integer, parameter :: NEX_XI_VAL = ',NEX_XI
+ write(IOUT,*) 'integer, parameter :: NEX_ETA_VAL = ',NEX_ETA
+ write(IOUT,*)
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE = ',NSPEC(IREGION_CRUST_MANTLE)
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE = ',NSPEC(IREGION_OUTER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE = ',NSPEC(IREGION_INNER_CORE)
+ write(IOUT,*)
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE = ',nglob(IREGION_CRUST_MANTLE)
+ write(IOUT,*) 'integer, parameter :: NGLOB_OUTER_CORE = ',nglob(IREGION_OUTER_CORE)
+ write(IOUT,*) 'integer, parameter :: NGLOB_INNER_CORE = ',nglob(IREGION_INNER_CORE)
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPECMAX_ANISO_IC = ',NSPECMAX_ANISO_IC
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPECMAX_ISO_MANTLE = ',NSPECMAX_ISO_MANTLE
+ write(IOUT,*) 'integer, parameter :: NSPECMAX_TISO_MANTLE = ',NSPECMAX_TISO_MANTLE
+ write(IOUT,*) 'integer, parameter :: NSPECMAX_ANISO_MANTLE = ',NSPECMAX_ANISO_MANTLE
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ATTENUAT = ',NSPEC_CRUST_MANTLE_ATTENUAT
+ write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_ATTENUATION = ',NSPEC_INNER_CORE_ATTENUATION
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STR_OR_ATT = ',NSPEC_CRUST_MANTLE_STR_OR_ATT
+ write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_STR_OR_ATT = ',NSPEC_INNER_CORE_STR_OR_ATT
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STR_AND_ATT = ',NSPEC_CRUST_MANTLE_STR_AND_ATT
+ write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_STR_AND_ATT = ',NSPEC_INNER_CORE_STR_AND_ATT
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STRAIN_ONLY = ',NSPEC_CRUST_MANTLE_STRAIN_ONLY
+ write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_STRAIN_ONLY = ',NSPEC_INNER_CORE_STRAIN_ONLY
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT = ',NSPEC_CRUST_MANTLE_ADJOINT
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ADJOINT = ',NSPEC_OUTER_CORE_ADJOINT
+ write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_ADJOINT = ',NSPEC_INNER_CORE_ADJOINT
+
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_ADJOINT = ',NGLOB_CRUST_MANTLE_ADJOINT
+ write(IOUT,*) 'integer, parameter :: NGLOB_OUTER_CORE_ADJOINT = ',NGLOB_OUTER_CORE_ADJOINT
+ write(IOUT,*) 'integer, parameter :: NGLOB_INNER_CORE_ADJOINT = ',NGLOB_INNER_CORE_ADJOINT
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ROT_ADJOINT = ',NSPEC_OUTER_CORE_ROT_ADJOINT
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STACEY = ',NSPEC_CRUST_MANTLE_STACEY
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_STACEY = ',NSPEC_OUTER_CORE_STACEY
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_OCEANS = ',NGLOB_CRUST_MANTLE_OCEANS
+ write(IOUT,*)
+
+! this to allow for code elimination by compiler in solver for performance
+
+ if(TRANSVERSE_ISOTROPY) then
+ write(IOUT,*) 'logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ANISOTROPIC_3D_MANTLE) then
+ write(IOUT,*) 'logical, parameter :: ANISOTROPIC_3D_MANTLE_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ANISOTROPIC_3D_MANTLE_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ANISOTROPIC_INNER_CORE) then
+ write(IOUT,*) 'logical, parameter :: ANISOTROPIC_INNER_CORE_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ANISOTROPIC_INNER_CORE_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ATTENUATION) then
+ write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ATTENUATION_NEW) then
+ write(IOUT,*) 'logical, parameter :: ATTENUATION_NEW_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ATTENUATION_NEW_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ATTENUATION_3D) then
+ write(IOUT,*) 'logical, parameter :: ATTENUATION_3D_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ATTENUATION_3D_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ELLIPTICITY) then
+ write(IOUT,*) 'logical, parameter :: ELLIPTICITY_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ELLIPTICITY_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(GRAVITY) then
+ write(IOUT,*) 'logical, parameter :: GRAVITY_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: GRAVITY_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(OCEANS) then
+ write(IOUT,*) 'logical, parameter :: OCEANS_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: OCEANS_VAL = .false.'
+ endif
+ write(IOUT,*)
+
+ if(ROTATION) then
+ write(IOUT,*) 'logical, parameter :: ROTATION_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: ROTATION_VAL = .false.'
+ endif
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ROTATION = ',NSPEC_OUTER_CORE_ROTATION
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NPROC_XI_VAL = ',NPROC_XI
+ write(IOUT,*) 'integer, parameter :: NPROC_ETA_VAL = ',NPROC_ETA
+ write(IOUT,*) 'integer, parameter :: NCHUNKS_VAL = ',NCHUNKS
+ write(IOUT,*) 'integer, parameter :: NPROCTOT_VAL = ',NPROCTOT
+ write(IOUT,*)
+
+ if(ATTENUATION) then
+ if(USE_3D_ATTENUATION_ARRAYS) then
+ att1 = NGLLX
+ att2 = NGLLY
+ att3 = NGLLZ
+ else
+ att1 = 1
+ att2 = 1
+ att3 = 1
+ endif
+ att4 = NSPEC(IREGION_CRUST_MANTLE)
+ att5 = NSPEC(IREGION_INNER_CORE)
+ else
+ att1 = 1
+ att2 = 1
+ att3 = 1
+ att4 = 1
+ att5 = 1
+ endif
+
+ write(IOUT,*) 'integer, parameter :: ATT1 = ',att1
+ write(IOUT,*) 'integer, parameter :: ATT2 = ',att2
+ write(IOUT,*) 'integer, parameter :: ATT3 = ',att3
+ write(IOUT,*) 'integer, parameter :: ATT4 = ',att4
+ write(IOUT,*) 'integer, parameter :: ATT5 = ',att5
+ write(IOUT,*)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_CM = ',NSPEC2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_CM = ',NSPEC2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_CM = ',NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_CM = ',NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_IC = ',NSPEC2DMAX_XMIN_XMAX(IREGION_INNER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_IC = ',NSPEC2DMAX_YMIN_YMAX(IREGION_INNER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_IC = ',NSPEC2D_BOTTOM(IREGION_INNER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_IC = ',NSPEC2D_TOP(IREGION_INNER_CORE)
+
+ write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_OC = ',NSPEC2DMAX_XMIN_XMAX(IREGION_OUTER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_OC = ',NSPEC2DMAX_YMIN_YMAX(IREGION_OUTER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_OC = ',NSPEC2D_BOTTOM(IREGION_OUTER_CORE)
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_OC = ',NSPEC2D_TOP(IREGION_OUTER_CORE)
+
+ ! for boundary kernels
+
+ if (SAVE_BOUNDARY_MESH) then
+ NSPEC2D_MOHO = NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ NSPEC2D_400 = NSPEC2D_MOHO / 4
+ NSPEC2D_670 = NSPEC2D_400
+ NSPEC2D_CMB = NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE)
+ NSPEC2D_ICB = NSPEC2D_BOTTOM(IREGION_OUTER_CORE)
+ else
+ NSPEC2D_MOHO = 1
+ NSPEC2D_400 = 1
+ NSPEC2D_670 = 1
+ NSPEC2D_CMB = 1
+ NSPEC2D_ICB = 1
+ endif
+
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_MOHO = ',NSPEC2D_MOHO
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_400 = ',NSPEC2D_400
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_670 = ',NSPEC2D_670
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_CMB = ',NSPEC2D_CMB
+ write(IOUT,*) 'integer, parameter :: NSPEC2D_ICB = ',NSPEC2D_ICB
+ write(IOUT,*)
+
+ ! deville routines only implemented for NGLLX = NGLLY = NGLLZ = 5
+ if( NGLLX == 5 .and. NGLLY == 5 .and. NGLLZ == 5 ) then
+ write(IOUT,*) 'logical, parameter :: USE_DEVILLE_PRODUCTS_VAL = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: USE_DEVILLE_PRODUCTS_VAL = .false.'
+ endif
+
+ ! backward/reconstruction of forward wavefield:
+ ! can only mimic attenuation effects on velocity at this point, since no full wavefield snapshots are stored
+ if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+
+ ! attenuation mimic:
+ ! mimicking effect of attenuation on apparent velocities, not amplitudes. that is,
+ ! phase shifts should be correctly accounted for, but amplitudes will differ in adjoint simulations
+
+!daniel: att - debug - check if mimic is still needed
+ if( ATTENUATION ) then
+ write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .true.'
+ else
+ write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
+ endif
+
+ else
+
+ ! calculates full attenuation (phase & amplitude effects) if used
+ write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
+ endif
+
+ ! attenuation and/or adjoint simulations
+ if (ATTENUATION .or. SIMULATION_TYPE /= 1 .or. SAVE_FORWARD &
+ .or. (MOVIE_VOLUME .and. SIMULATION_TYPE /= 3)) then
+ write(IOUT,*) 'logical, parameter :: COMPUTE_AND_STORE_STRAIN = .true. '
+ else
+ write(IOUT,*) 'logical, parameter :: COMPUTE_AND_STORE_STRAIN = .false.'
+ endif
+ write(IOUT,*)
+
+ if (MOVIE_VOLUME) then
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = NSPEC_CRUST_MANTLE'
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = NGLOB_CRUST_MANTLE'
+ else
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = 1'
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = 1'
+ endif
+
+ close(IOUT)
+
+ end subroutine save_header_file
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -1,544 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Princeton University, USA
-! and CNRS / INRIA / University of Pau, France
-! (c) Princeton University and CNRS / INRIA / University of Pau
-! April 2011
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-! save header file OUTPUT_FILES/values_from_mesher.h
-
- subroutine save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
- TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION, &
- OCEANS,ATTENUATION,ATTENUATION_NEW,ATTENUATION_3D, &
- ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
- INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES, &
- CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP,&
- static_memory_size, &
- NSPEC2D_TOP,NSPEC2D_BOTTOM, &
- NSPEC2DMAX_YMIN_YMAX,NSPEC2DMAX_XMIN_XMAX, &
- NPROC_XI,NPROC_ETA, &
- NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
- NSPEC_INNER_CORE_ATTENUATION, &
- NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
- NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
- NSPEC_CRUST_MANTLE_ADJOINT, &
- NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
- NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
- NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
- NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
- SIMULATION_TYPE,SAVE_FORWARD,MOVIE_VOLUME)
-
- implicit none
-
- include "constants.h"
-
- integer, dimension(MAX_NUM_REGIONS) :: NSPEC, nglob
-
- integer NEX_XI,NEX_ETA,NPROC,NPROCTOT,NCHUNKS,NSOURCES,NSTEP
-
- logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_NEW,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE
-
- double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
- CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH
-
- ! static memory size needed by the solver
- double precision :: static_memory_size
-
- integer, dimension(MAX_NUM_REGIONS) :: &
- NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_YMIN_YMAX,NSPEC2DMAX_XMIN_XMAX
- integer :: NPROC_XI,NPROC_ETA
-
- integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
- NSPEC_INNER_CORE_ATTENUATION, &
- NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
- NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
- NSPEC_CRUST_MANTLE_ADJOINT, &
- NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
- NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
- NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
- NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
- NSPEC2D_MOHO, NSPEC2D_400, NSPEC2D_670, NSPEC2D_CMB, NSPEC2D_ICB
-
- integer :: SIMULATION_TYPE
- logical :: SAVE_FORWARD,MOVIE_VOLUME
-
- ! local parameters
- double precision :: subtract_central_cube_elems,subtract_central_cube_points
- ! for regional code
- double precision x,y,gamma,rgt,xi,eta
- double precision x_top,y_top,z_top
- double precision ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD
- ! rotation matrix from Euler angles
- integer i,j,ix,iy,icorner
- double precision rotation_matrix(3,3)
- double precision vector_ori(3),vector_rotated(3)
- double precision r_corner,theta_corner,phi_corner,lat,long,colat_corner
- integer :: att1,att2,att3,att4,att5
- integer :: ier
- character(len=150) HEADER_FILE
-
- ! copy number of elements and points in an include file for the solver
- call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
- open(unit=IOUT,file=HEADER_FILE,status='unknown',iostat=ier)
- if( ier /= 0 ) stop 'error opening OUTPUT_FILES/values_from_mesher.h'
-
- write(IOUT,*)
-
- write(IOUT,*) '!'
- write(IOUT,*) '! this is the parameter file for static compilation of the solver'
- write(IOUT,*) '!'
- write(IOUT,*) '! mesh statistics:'
- write(IOUT,*) '! ---------------'
- write(IOUT,*) '!'
- write(IOUT,*) '!'
- write(IOUT,*) '! number of chunks = ',NCHUNKS
- write(IOUT,*) '!'
-
-! the central cube is counted 6 times, therefore remove 5 times
- if(INCLUDE_CENTRAL_CUBE) then
- write(IOUT,*) '! these statistics include the central cube'
- subtract_central_cube_elems = 5.d0 * dble((NEX_XI/8))**3
- subtract_central_cube_points = 5.d0 * (dble(NEX_XI/8)*dble(NGLLX-1)+1.d0)**3
- else
- write(IOUT,*) '! these statistics do not include the central cube'
- subtract_central_cube_elems = 0.d0
- subtract_central_cube_points = 0.d0
- endif
-
- write(IOUT,*) '!'
- write(IOUT,*) '! number of processors = ',NPROCTOT ! should be = NPROC
- write(IOUT,*) '!'
- write(IOUT,*) '! maximum number of points per region = ',nglob(IREGION_CRUST_MANTLE)
- write(IOUT,*) '!'
-! use fused loops on NEC SX
- write(IOUT,*) '! on NEC SX, make sure "loopcnt=" parameter'
- write(IOUT,*) '! in Makefile is greater than max vector length = ',nglob(IREGION_CRUST_MANTLE)*NDIM
- write(IOUT,*) '!'
-
- write(IOUT,*) '! total elements per slice = ',sum(NSPEC)
- write(IOUT,*) '! total points per slice = ',sum(nglob)
- write(IOUT,*) '!'
-
- write(IOUT,'(1x,a,i1,a)') '! total for full ',NCHUNKS,'-chunk mesh:'
- write(IOUT,*) '! ---------------------------'
- write(IOUT,*) '!'
- write(IOUT,*) '! exact total number of spectral elements in entire mesh = '
- write(IOUT,*) '! ',dble(NCHUNKS)*dble(NPROC)*dble(sum(NSPEC)) - subtract_central_cube_elems
- write(IOUT,*) '! approximate total number of points in entire mesh = '
- write(IOUT,*) '! ',dble(NCHUNKS)*dble(NPROC)*dble(sum(nglob)) - subtract_central_cube_points
-! there are 3 DOFs in solid regions, but only 1 in fluid outer core
- write(IOUT,*) '! approximate total number of degrees of freedom in entire mesh = '
- write(IOUT,*) '! ',dble(NCHUNKS)*dble(NPROC)*(3.d0*(dble(sum(nglob))) &
- - 2.d0*dble(nglob(IREGION_OUTER_CORE))) &
- - 3.d0*subtract_central_cube_points
- write(IOUT,*) '!'
-
-! display location of chunk if regional run
- if(NCHUNKS /= 6) then
-
- write(IOUT,*) '! position of the mesh chunk at the surface:'
- write(IOUT,*) '! -----------------------------------------'
- write(IOUT,*) '!'
- write(IOUT,*) '! angular size in first direction in degrees = ',sngl(ANGULAR_WIDTH_XI_IN_DEGREES)
- write(IOUT,*) '! angular size in second direction in degrees = ',sngl(ANGULAR_WIDTH_ETA_IN_DEGREES)
- write(IOUT,*) '!'
- write(IOUT,*) '! longitude of center in degrees = ',sngl(CENTER_LONGITUDE_IN_DEGREES)
- write(IOUT,*) '! latitude of center in degrees = ',sngl(CENTER_LATITUDE_IN_DEGREES)
- write(IOUT,*) '!'
- write(IOUT,*) '! angle of rotation of the first chunk = ',sngl(GAMMA_ROTATION_AZIMUTH)
-
-! convert width to radians
- ANGULAR_WIDTH_XI_RAD = ANGULAR_WIDTH_XI_IN_DEGREES * DEGREES_TO_RADIANS
- ANGULAR_WIDTH_ETA_RAD = ANGULAR_WIDTH_ETA_IN_DEGREES * DEGREES_TO_RADIANS
-
-! compute rotation matrix from Euler angles
- call euler_angles(rotation_matrix,CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH)
-
-! loop on the four corners of the chunk to display their coordinates
- icorner = 0
- do iy = 0,1
- do ix = 0,1
-
- icorner = icorner + 1
-
- xi= - ANGULAR_WIDTH_XI_RAD/2. + dble(ix)*ANGULAR_WIDTH_XI_RAD
- eta= - ANGULAR_WIDTH_ETA_RAD/2. + dble(iy)*ANGULAR_WIDTH_ETA_RAD
-
- x=dtan(xi)
- y=dtan(eta)
-
- gamma=ONE/dsqrt(ONE+x*x+y*y)
- rgt=R_UNIT_SPHERE*gamma
-
-! define the mesh points at the top surface
- x_top=-y*rgt
- y_top=x*rgt
- z_top=rgt
-
-! rotate top
- vector_ori(1) = x_top
- vector_ori(2) = y_top
- vector_ori(3) = z_top
- do i=1,3
- vector_rotated(i)=0.0d0
- do j=1,3
- vector_rotated(i)=vector_rotated(i)+rotation_matrix(i,j)*vector_ori(j)
- enddo
- enddo
- x_top = vector_rotated(1)
- y_top = vector_rotated(2)
- z_top = vector_rotated(3)
-
-! convert to latitude and longitude
- call xyz_2_rthetaphi_dble(x_top,y_top,z_top,r_corner,theta_corner,phi_corner)
- call reduce(theta_corner,phi_corner)
-
-! convert geocentric to geographic colatitude
- colat_corner=PI_OVER_TWO-datan(1.006760466d0*dcos(theta_corner)/dmax1(TINYVAL,dsin(theta_corner)))
- if(phi_corner>PI) phi_corner=phi_corner-TWO_PI
-
-! compute real position of the source
- lat = (PI_OVER_TWO-colat_corner)*RADIANS_TO_DEGREES
- long = phi_corner*RADIANS_TO_DEGREES
-
- write(IOUT,*) '!'
- write(IOUT,*) '! corner ',icorner
- write(IOUT,*) '! longitude in degrees = ',long
- write(IOUT,*) '! latitude in degrees = ',lat
-
- enddo
- enddo
-
- write(IOUT,*) '!'
-
- endif ! regional chunk
-
- write(IOUT,*) '! resolution of the mesh at the surface:'
- write(IOUT,*) '! -------------------------------------'
- write(IOUT,*) '!'
- write(IOUT,*) '! spectral elements along a great circle = ',4*NEX_XI
- write(IOUT,*) '! GLL points along a great circle = ',4*NEX_XI*(NGLLX-1)
- write(IOUT,*) '! average distance between points in degrees = ',360./real(4*NEX_XI*(NGLLX-1))
- write(IOUT,*) '! average distance between points in km = ',real(TWO_PI*R_EARTH/1000.d0)/real(4*NEX_XI*(NGLLX-1))
- write(IOUT,*) '! average size of a spectral element in km = ',real(TWO_PI*R_EARTH/1000.d0)/real(4*NEX_XI)
- write(IOUT,*) '!'
- write(IOUT,*) '! number of time steps = ',NSTEP
- write(IOUT,*) '!'
- write(IOUT,*) '! number of seismic sources = ',NSOURCES
- write(IOUT,*) '!'
- write(IOUT,*)
-
- write(IOUT,*) '! approximate static memory needed by the solver:'
- write(IOUT,*) '! ----------------------------------------------'
- write(IOUT,*) '!'
- write(IOUT,*) '! (lower bound, usually the real amount used is 5% to 10% higher)'
- write(IOUT,*) '!'
- write(IOUT,*) '! (you can get a more precise estimate of the size used per MPI process'
- write(IOUT,*) '! by typing "size -d bin/xspecfem3D"'
- write(IOUT,*) '! after compiling the code with the DATA/Par_file you plan to use)'
- write(IOUT,*) '!'
- write(IOUT,*) '! size of static arrays per slice = ',static_memory_size/1.d6,' MB'
- write(IOUT,*) '! = ',static_memory_size/1048576.d0,' MiB'
- write(IOUT,*) '! = ',static_memory_size/1.d9,' GB'
- write(IOUT,*) '! = ',static_memory_size/1073741824.d0,' GiB'
- write(IOUT,*) '!'
-
- ! note: using less memory becomes an issue only if the strong scaling of the code is poor.
- ! Some users will run simulations with an executable using far less than 80% RAM per core
- ! if they prefer having a faster computational time (and use a higher number of cores).
-
- write(IOUT,*) '! (should be below to 80% or 90% of the memory installed per core)'
- write(IOUT,*) '! (if significantly more, the job will not run by lack of memory )'
- write(IOUT,*) '! (note that if significantly less, you waste a significant amount'
- write(IOUT,*) '! of memory per processor core)'
- write(IOUT,*) '! (but that can be perfectly acceptable if you can afford it and'
- write(IOUT,*) '! want faster results by using more cores)'
- write(IOUT,*) '!'
- if(static_memory_size*dble(NPROCTOT)/1.d6 < 10000.d0) then
- write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1.d6,' MB'
- write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1048576.d0,' MiB'
- write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1.d9,' GB'
- else
- write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1.d9,' GB'
- endif
- write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GiB'
- write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1.d12,' TB'
- write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TiB'
- write(IOUT,*) '!'
-
- write(IOUT,*)
- write(IOUT,*) 'integer, parameter :: NEX_XI_VAL = ',NEX_XI
- write(IOUT,*) 'integer, parameter :: NEX_ETA_VAL = ',NEX_ETA
- write(IOUT,*)
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE = ',NSPEC(IREGION_CRUST_MANTLE)
- write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE = ',NSPEC(IREGION_OUTER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE = ',NSPEC(IREGION_INNER_CORE)
- write(IOUT,*)
- write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE = ',nglob(IREGION_CRUST_MANTLE)
- write(IOUT,*) 'integer, parameter :: NGLOB_OUTER_CORE = ',nglob(IREGION_OUTER_CORE)
- write(IOUT,*) 'integer, parameter :: NGLOB_INNER_CORE = ',nglob(IREGION_INNER_CORE)
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPECMAX_ANISO_IC = ',NSPECMAX_ANISO_IC
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPECMAX_ISO_MANTLE = ',NSPECMAX_ISO_MANTLE
- write(IOUT,*) 'integer, parameter :: NSPECMAX_TISO_MANTLE = ',NSPECMAX_TISO_MANTLE
- write(IOUT,*) 'integer, parameter :: NSPECMAX_ANISO_MANTLE = ',NSPECMAX_ANISO_MANTLE
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ATTENUAT = ',NSPEC_CRUST_MANTLE_ATTENUAT
- write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_ATTENUATION = ',NSPEC_INNER_CORE_ATTENUATION
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STR_OR_ATT = ',NSPEC_CRUST_MANTLE_STR_OR_ATT
- write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_STR_OR_ATT = ',NSPEC_INNER_CORE_STR_OR_ATT
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STR_AND_ATT = ',NSPEC_CRUST_MANTLE_STR_AND_ATT
- write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_STR_AND_ATT = ',NSPEC_INNER_CORE_STR_AND_ATT
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STRAIN_ONLY = ',NSPEC_CRUST_MANTLE_STRAIN_ONLY
- write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_STRAIN_ONLY = ',NSPEC_INNER_CORE_STRAIN_ONLY
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT = ',NSPEC_CRUST_MANTLE_ADJOINT
- write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ADJOINT = ',NSPEC_OUTER_CORE_ADJOINT
- write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_ADJOINT = ',NSPEC_INNER_CORE_ADJOINT
-
- write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_ADJOINT = ',NGLOB_CRUST_MANTLE_ADJOINT
- write(IOUT,*) 'integer, parameter :: NGLOB_OUTER_CORE_ADJOINT = ',NGLOB_OUTER_CORE_ADJOINT
- write(IOUT,*) 'integer, parameter :: NGLOB_INNER_CORE_ADJOINT = ',NGLOB_INNER_CORE_ADJOINT
-
- write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ROT_ADJOINT = ',NSPEC_OUTER_CORE_ROT_ADJOINT
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_STACEY = ',NSPEC_CRUST_MANTLE_STACEY
- write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_STACEY = ',NSPEC_OUTER_CORE_STACEY
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_OCEANS = ',NGLOB_CRUST_MANTLE_OCEANS
- write(IOUT,*)
-
-! this to allow for code elimination by compiler in solver for performance
-
- if(TRANSVERSE_ISOTROPY) then
- write(IOUT,*) 'logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ANISOTROPIC_3D_MANTLE) then
- write(IOUT,*) 'logical, parameter :: ANISOTROPIC_3D_MANTLE_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ANISOTROPIC_3D_MANTLE_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ANISOTROPIC_INNER_CORE) then
- write(IOUT,*) 'logical, parameter :: ANISOTROPIC_INNER_CORE_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ANISOTROPIC_INNER_CORE_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ATTENUATION) then
- write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ATTENUATION_NEW) then
- write(IOUT,*) 'logical, parameter :: ATTENUATION_NEW_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ATTENUATION_NEW_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ATTENUATION_3D) then
- write(IOUT,*) 'logical, parameter :: ATTENUATION_3D_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ATTENUATION_3D_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ELLIPTICITY) then
- write(IOUT,*) 'logical, parameter :: ELLIPTICITY_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ELLIPTICITY_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(GRAVITY) then
- write(IOUT,*) 'logical, parameter :: GRAVITY_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: GRAVITY_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(OCEANS) then
- write(IOUT,*) 'logical, parameter :: OCEANS_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: OCEANS_VAL = .false.'
- endif
- write(IOUT,*)
-
- if(ROTATION) then
- write(IOUT,*) 'logical, parameter :: ROTATION_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: ROTATION_VAL = .false.'
- endif
- write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ROTATION = ',NSPEC_OUTER_CORE_ROTATION
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NPROC_XI_VAL = ',NPROC_XI
- write(IOUT,*) 'integer, parameter :: NPROC_ETA_VAL = ',NPROC_ETA
- write(IOUT,*) 'integer, parameter :: NCHUNKS_VAL = ',NCHUNKS
- write(IOUT,*) 'integer, parameter :: NPROCTOT_VAL = ',NPROCTOT
- write(IOUT,*)
-
- if(ATTENUATION) then
- if(USE_3D_ATTENUATION_ARRAYS) then
- att1 = NGLLX
- att2 = NGLLY
- att3 = NGLLZ
- else
- att1 = 1
- att2 = 1
- att3 = 1
- endif
- att4 = NSPEC(IREGION_CRUST_MANTLE)
- att5 = NSPEC(IREGION_INNER_CORE)
- else
- att1 = 1
- att2 = 1
- att3 = 1
- att4 = 1
- att5 = 1
- endif
-
- write(IOUT,*) 'integer, parameter :: ATT1 = ',att1
- write(IOUT,*) 'integer, parameter :: ATT2 = ',att2
- write(IOUT,*) 'integer, parameter :: ATT3 = ',att3
- write(IOUT,*) 'integer, parameter :: ATT4 = ',att4
- write(IOUT,*) 'integer, parameter :: ATT5 = ',att5
- write(IOUT,*)
-
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_CM = ',NSPEC2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE)
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_CM = ',NSPEC2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE)
- write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_CM = ',NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE)
- write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_CM = ',NSPEC2D_TOP(IREGION_CRUST_MANTLE)
-
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_IC = ',NSPEC2DMAX_XMIN_XMAX(IREGION_INNER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_IC = ',NSPEC2DMAX_YMIN_YMAX(IREGION_INNER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_IC = ',NSPEC2D_BOTTOM(IREGION_INNER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_IC = ',NSPEC2D_TOP(IREGION_INNER_CORE)
-
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_OC = ',NSPEC2DMAX_XMIN_XMAX(IREGION_OUTER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_OC = ',NSPEC2DMAX_YMIN_YMAX(IREGION_OUTER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_OC = ',NSPEC2D_BOTTOM(IREGION_OUTER_CORE)
- write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_OC = ',NSPEC2D_TOP(IREGION_OUTER_CORE)
-
- ! for boundary kernels
-
- if (SAVE_BOUNDARY_MESH) then
- NSPEC2D_MOHO = NSPEC2D_TOP(IREGION_CRUST_MANTLE)
- NSPEC2D_400 = NSPEC2D_MOHO / 4
- NSPEC2D_670 = NSPEC2D_400
- NSPEC2D_CMB = NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE)
- NSPEC2D_ICB = NSPEC2D_BOTTOM(IREGION_OUTER_CORE)
- else
- NSPEC2D_MOHO = 1
- NSPEC2D_400 = 1
- NSPEC2D_670 = 1
- NSPEC2D_CMB = 1
- NSPEC2D_ICB = 1
- endif
-
- write(IOUT,*) 'integer, parameter :: NSPEC2D_MOHO = ',NSPEC2D_MOHO
- write(IOUT,*) 'integer, parameter :: NSPEC2D_400 = ',NSPEC2D_400
- write(IOUT,*) 'integer, parameter :: NSPEC2D_670 = ',NSPEC2D_670
- write(IOUT,*) 'integer, parameter :: NSPEC2D_CMB = ',NSPEC2D_CMB
- write(IOUT,*) 'integer, parameter :: NSPEC2D_ICB = ',NSPEC2D_ICB
- write(IOUT,*)
-
- ! deville routines only implemented for NGLLX = NGLLY = NGLLZ = 5
- if( NGLLX == 5 .and. NGLLY == 5 .and. NGLLZ == 5 ) then
- write(IOUT,*) 'logical, parameter :: USE_DEVILLE_PRODUCTS_VAL = .true.'
- else
- write(IOUT,*) 'logical, parameter :: USE_DEVILLE_PRODUCTS_VAL = .false.'
- endif
-
- ! backward/reconstruction of forward wavefield:
- ! can only mimic attenuation effects on velocity at this point, since no full wavefield snapshots are stored
- if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
-
- ! attenuation mimic:
- ! mimicking effect of attenuation on apparent velocities, not amplitudes. that is,
- ! phase shifts should be correctly accounted for, but amplitudes will differ in adjoint simulations
-
-!daniel: att - debug - check if mimic is still needed
- if( ATTENUATION ) then
- write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .true.'
- else
- write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
- endif
-
- else
-
- ! calculates full attenuation (phase & amplitude effects) if used
- write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
- endif
-
- ! attenuation and/or adjoint simulations
- if (ATTENUATION .or. SIMULATION_TYPE /= 1 .or. SAVE_FORWARD &
- .or. (MOVIE_VOLUME .and. SIMULATION_TYPE /= 3)) then
- write(IOUT,*) 'logical, parameter :: COMPUTE_AND_STORE_STRAIN = .true. '
- else
- write(IOUT,*) 'logical, parameter :: COMPUTE_AND_STORE_STRAIN = .false.'
- endif
- write(IOUT,*)
-
- if (MOVIE_VOLUME) then
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = NSPEC_CRUST_MANTLE'
- write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = NGLOB_CRUST_MANTLE'
- else
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = 1'
- write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = 1'
- endif
-
- close(IOUT)
-
- end subroutine save_header_file
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_classical.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_classical.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -7,6 +7,15 @@
! Newmark time scheme update
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ displ_crust_mantle(i,1) = displ_crust_mantle(i,1) &
+ + deltat*veloc_crust_mantle(i,1) + deltatsqover2*accel_crust_mantle(i,1)
+ veloc_crust_mantle(i,1) = veloc_crust_mantle(i,1) &
+ + deltatover2*accel_crust_mantle(i,1)
+ accel_crust_mantle(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
displ_crust_mantle(:,i) = displ_crust_mantle(:,i) &
+ deltat*veloc_crust_mantle(:,i) + deltatsqover2*accel_crust_mantle(:,i)
@@ -14,6 +23,8 @@
+ deltatover2*accel_crust_mantle(:,i)
accel_crust_mantle(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
! outer core
do i=1,NGLOB_OUTER_CORE
displ_outer_core(i) = displ_outer_core(i) &
@@ -22,7 +33,17 @@
+ deltatover2*accel_outer_core(i)
accel_outer_core(i) = 0._CUSTOM_REAL
enddo
+
! inner core
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ displ_inner_core(i,1) = displ_inner_core(i,1) &
+ + deltat*veloc_inner_core(i,1) + deltatsqover2*accel_inner_core(i,1)
+ veloc_inner_core(i,1) = veloc_inner_core(i,1) &
+ + deltatover2*accel_inner_core(i,1)
+ accel_inner_core(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_INNER_CORE
displ_inner_core(:,i) = displ_inner_core(:,i) &
+ deltat*veloc_inner_core(:,i) + deltatsqover2*accel_inner_core(:,i)
@@ -30,6 +51,7 @@
+ deltatover2*accel_inner_core(:,i)
accel_inner_core(:,i) = 0._CUSTOM_REAL
enddo
+ endif
! integral of strain for adjoint movie volume
if(MOVIE_VOLUME .and. (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) ) then
@@ -92,9 +114,10 @@
buffer_send_faces,buffer_received_faces,npoin2D_max_all_CM_IC, &
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -117,7 +140,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
! Stacey absorbing boundaries
@@ -218,9 +241,10 @@
buffer_send_faces,buffer_received_faces,npoin2D_max_all_CM_IC, &
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -242,7 +266,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
do while (iphase <= 7) ! make sure the last communications are finished and processed
@@ -267,7 +291,7 @@
veloc_outer_core(i) = veloc_outer_core(i) + deltatover2*accel_outer_core(i)
enddo
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
! ****************************************************
! big loop over all spectral elements in the solid
@@ -306,7 +330,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -317,11 +341,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -358,11 +382,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
! Deville routine
@@ -390,16 +414,16 @@
npoin2D_cube_from_slices,buffer_all_cube_from_slices,buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -428,17 +452,16 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
! Stacey
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
call compute_stacey_crust_mantle_forward(ichunk, &
it,SAVE_FORWARD,ibool_crust_mantle, &
veloc_crust_mantle,accel_crust_mantle, &
@@ -469,7 +492,7 @@
accel_crust_mantle,sourcearrays, &
DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
islice_selected_source,ispec_selected_source,it, &
- hdur,xi_source,eta_source,gamma_source,nu_source,istage)
+ hdur,xi_source,eta_source,gamma_source,nu_source,istage,USE_LDDRK)
! add adjoint sources only if adjoint simulation is performed for source inversion only
!! DK DK UNDO_ATTENUATION this must remain here even when SIMULATION_TYPE == 3 because it applies to array
@@ -486,14 +509,6 @@
it,it_begin,station_name,network_name,DT)
endif
-! ! add adjoint sources and add sources for backward/reconstructed wavefield
-! if (SIMULATION_TYPE == 3) &
-! call compute_add_sources_backward(myrank,NSOURCES,NSTEP, &
-! b_accel_crust_mantle,sourcearrays, &
-! DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
-! islice_selected_source,ispec_selected_source,it, &
-! hdur,xi_source,eta_source,gamma_source,nu_source)
-
! NOISE_TOMOGRAPHY
if ( NOISE_TOMOGRAPHY == 1 ) then
! the first step of noise tomography is to use |S(\omega)|^2 as a point force source at one of the receivers.
@@ -610,7 +625,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -621,11 +636,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -662,11 +677,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
! Deville routine
@@ -694,16 +709,16 @@
npoin2D_cube_from_slices,buffer_all_cube_from_slices,buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -732,12 +747,12 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
! assemble all the contributions between slices using MPI
@@ -775,7 +790,6 @@
endif ! end of assembling forces with the central cube
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
do i=1,NGLOB_CRUST_MANTLE
accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
@@ -783,9 +797,7 @@
- two_omega_earth*veloc_crust_mantle(1,i)
accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
enddo
-
else
-
do i=1,NGLOB_CRUST_MANTLE
accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
@@ -793,10 +805,9 @@
- two_omega_earth*veloc_crust_mantle(1,i)
accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
enddo
-
endif
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
! couples ocean with crust mantle
if(OCEANS_VAL) then
@@ -804,9 +815,9 @@
rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
endif
!-------------------------------------------------------------------------------------------------
@@ -819,9 +830,16 @@
! Newmark time scheme - corrector for elastic parts
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ veloc_crust_mantle(i,1) = veloc_crust_mantle(i,1) + deltatover2*accel_crust_mantle(i,1)
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
enddo
+ endif
+
! inner core
do i=1,NGLOB_INNER_CORE
accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
@@ -829,9 +847,17 @@
accel_inner_core(2,i) = accel_inner_core(2,i)*rmass_inner_core(i) &
- two_omega_earth*veloc_inner_core(1,i)
accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ enddo
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ veloc_inner_core(i,1) = veloc_inner_core(i,1) + deltatover2*accel_inner_core(i,1)
+ enddo
+ else
+ do i=1,NGLOB_INNER_CORE
veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
enddo
+ endif
! write the seismograms with time shift
@@ -948,7 +974,8 @@
mask_3dmovie,nu_3dmovie)
else if (MOVIE_VOLUME_TYPE == 4) then ! output divergence and curl in whole volume
-!!!!! for undo_att this type of MOVIE is not supported
+ stop 'MOVIE_VOLUME_TYPE == 4 is not implemented for UNDO_ATTENUATION in order to save memory'
+!!!!! for UNDO_ATTENUATION this type of MOVIE is not supported
!!! call write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
!!! div_displ_outer_core, &
!!! accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
@@ -982,4 +1009,3 @@
endif
endif ! MOVIE_VOLUME
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_undo_att.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part1_undo_att.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -4,19 +4,31 @@
!! DK DK it also handles the cases NOISE_TOMOGRAPHY == 1 and NOISE_TOMOGRAPHY == 2
!! DK DK
- ! Newmark time scheme update
+ ! Newmark or LDDRK time scheme update
+ do istage = 1, NSTAGE_TIME_SCHEME ! is equal to 1 if Newmark because only one stage then
- do istage = 1, NSTAGE_TIME_SCHEME !ZN begin loop of istage
if(USE_LDDRK)then
+
! mantle
accel_crust_mantle(:,:) = 0._CUSTOM_REAL
! outer core
accel_outer_core(:) = 0._CUSTOM_REAL
! inner core
accel_inner_core(:,:) = 0._CUSTOM_REAL
+
else
+
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ displ_crust_mantle(i,1) = displ_crust_mantle(i,1) &
+ + deltat*veloc_crust_mantle(i,1) + deltatsqover2*accel_crust_mantle(i,1)
+ veloc_crust_mantle(i,1) = veloc_crust_mantle(i,1) &
+ + deltatover2*accel_crust_mantle(i,1)
+ accel_crust_mantle(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
displ_crust_mantle(:,i) = displ_crust_mantle(:,i) &
+ deltat*veloc_crust_mantle(:,i) + deltatsqover2*accel_crust_mantle(:,i)
@@ -24,6 +36,8 @@
+ deltatover2*accel_crust_mantle(:,i)
accel_crust_mantle(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
! outer core
do i=1,NGLOB_OUTER_CORE
displ_outer_core(i) = displ_outer_core(i) &
@@ -32,7 +46,17 @@
+ deltatover2*accel_outer_core(i)
accel_outer_core(i) = 0._CUSTOM_REAL
enddo
+
! inner core
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ displ_inner_core(i,1) = displ_inner_core(i,1) &
+ + deltat*veloc_inner_core(i,1) + deltatsqover2*accel_inner_core(i,1)
+ veloc_inner_core(i,1) = veloc_inner_core(i,1) &
+ + deltatover2*accel_inner_core(i,1)
+ accel_inner_core(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_INNER_CORE
displ_inner_core(:,i) = displ_inner_core(:,i) &
+ deltat*veloc_inner_core(:,i) + deltatsqover2*accel_inner_core(:,i)
@@ -40,6 +64,8 @@
+ deltatover2*accel_inner_core(:,i)
accel_inner_core(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
endif
if(istage == 1)then
@@ -112,9 +138,10 @@
buffer_send_faces,buffer_received_faces,npoin2D_max_all_CM_IC, &
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -137,7 +164,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
! Stacey absorbing boundaries
@@ -238,9 +265,10 @@
buffer_send_faces,buffer_received_faces,npoin2D_max_all_CM_IC, &
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -262,7 +290,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
do while (iphase <= 7) ! make sure the last communications are finished and processed
@@ -297,7 +325,8 @@
enddo
endif
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
+
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -335,7 +364,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -346,11 +375,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -387,11 +416,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
! Deville routine
@@ -419,16 +448,16 @@
npoin2D_cube_from_slices,buffer_all_cube_from_slices,buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -457,17 +486,16 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
! Stacey
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
call compute_stacey_crust_mantle_forward(ichunk, &
it,SAVE_FORWARD,ibool_crust_mantle, &
veloc_crust_mantle,accel_crust_mantle, &
@@ -498,7 +526,7 @@
accel_crust_mantle,sourcearrays, &
DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
islice_selected_source,ispec_selected_source,it, &
- hdur,xi_source,eta_source,gamma_source,nu_source,istage)
+ hdur,xi_source,eta_source,gamma_source,nu_source,istage,USE_LDDRK)
! add adjoint sources only if adjoint simulation is performed for source inversion only
!! DK DK UNDO_ATTENUATION this must remain here even when SIMULATION_TYPE == 3 because it applies to array
@@ -627,7 +655,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -638,11 +666,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -679,11 +707,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
! Deville routine
@@ -711,16 +739,16 @@
npoin2D_cube_from_slices,buffer_all_cube_from_slices,buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -749,12 +777,12 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
! assemble all the contributions between slices using MPI
@@ -791,8 +819,8 @@
enddo
endif ! end of assembling forces with the central cube
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
+ if((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) .and. .not. USE_LDDRK) then
do i=1,NGLOB_CRUST_MANTLE
accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
@@ -800,9 +828,7 @@
- two_omega_earth*veloc_crust_mantle(1,i)
accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
enddo
-
else
-
do i=1,NGLOB_CRUST_MANTLE
accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
@@ -810,19 +836,19 @@
- two_omega_earth*veloc_crust_mantle(1,i)
accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
enddo
-
endif
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
+
! couples ocean with crust mantle
if(OCEANS_VAL) then
call compute_coupling_ocean(accel_crust_mantle, &
rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
endif
!-------------------------------------------------------------------------------------------------
@@ -835,7 +861,7 @@
! Newmark time scheme - corrector for elastic parts
! mantle
- if(USE_LDDRK)then
+ if(USE_LDDRK) then
do i=1,NGLOB_CRUST_MANTLE
veloc_crust_mantle_lddrk(:,i) = ALPHA_LDDRK(istage) * veloc_crust_mantle_lddrk(:,i) &
+ deltat * accel_crust_mantle(:,i)
@@ -845,31 +871,59 @@
displ_crust_mantle(:,i) = displ_crust_mantle(:,i) + BETA_LDDRK(istage) * displ_crust_mantle_lddrk(:,i)
enddo
else
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ veloc_crust_mantle(i,1) = veloc_crust_mantle(i,1) + deltatover2*accel_crust_mantle(i,1)
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
enddo
+ endif
endif
+
! inner core
- do i=1,NGLOB_INNER_CORE
- accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
- + two_omega_earth*veloc_inner_core(2,i)
- accel_inner_core(2,i) = accel_inner_core(2,i)*rmass_inner_core(i) &
- - two_omega_earth*veloc_inner_core(1,i)
- accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. .not. USE_LDDRK) then
+ do i=1,NGLOB_INNER_CORE
+ accel_inner_core(1,i) = accel_inner_core(1,i)*rmassx_inner_core(i) &
+ + two_omega_earth*veloc_inner_core(2,i)
+ accel_inner_core(2,i) = accel_inner_core(2,i)*rmassy_inner_core(i) &
+ - two_omega_earth*veloc_inner_core(1,i)
+ accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ enddo
+ else
+ do i=1,NGLOB_INNER_CORE
+ accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
+ + two_omega_earth*veloc_inner_core(2,i)
+ accel_inner_core(2,i) = accel_inner_core(2,i)*rmass_inner_core(i) &
+ - two_omega_earth*veloc_inner_core(1,i)
+ accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ enddo
+ endif
- if(USE_LDDRK)then
+ if(USE_LDDRK)then
+ do i=1,NGLOB_INNER_CORE
veloc_inner_core_lddrk(:,i) = ALPHA_LDDRK(istage) * veloc_inner_core_lddrk(:,i) &
+ deltat * accel_inner_core(:,i)
displ_inner_core_lddrk(:,i) = ALPHA_LDDRK(istage) * displ_inner_core_lddrk(:,i) &
+ deltat * veloc_inner_core(:,i)
veloc_inner_core(:,i) = veloc_inner_core(:,i) + BETA_LDDRK(istage) * veloc_inner_core_lddrk(:,i)
displ_inner_core(:,i) = displ_inner_core(:,i) + BETA_LDDRK(istage) * displ_inner_core_lddrk(:,i)
- else
+ enddo
+ else
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ veloc_inner_core(i,1) = veloc_inner_core(i,1) + deltatover2*accel_inner_core(i,1)
+ enddo
+ else
+ do i=1,NGLOB_INNER_CORE
veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
- endif
- enddo
- enddo !ZN end loop of istage
+ enddo
+ endif
+ endif
+ enddo ! end of very big external loop on istage for all the stages of the LDDRK time scheme (only one stage if Newmark)
+
! write the seismograms with time shift
! store the seismograms only if there is at least one receiver located in this slice
@@ -1015,7 +1069,8 @@
mask_3dmovie,nu_3dmovie)
else if (MOVIE_VOLUME_TYPE == 4) then ! output divergence and curl in whole volume
-!!!!!! for undo_att this type of MOVIE is not supported
+ stop 'MOVIE_VOLUME_TYPE == 4 is not implemented for UNDO_ATTENUATION in order to save memory'
+!!!!! for UNDO_ATTENUATION this type of MOVIE is not supported
!!! call write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
!!! div_displ_outer_core, &
!!! accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
@@ -1049,4 +1104,3 @@
endif
endif ! MOVIE_VOLUME
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_classical.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_classical.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -8,7 +8,17 @@
! backward field
if (SIMULATION_TYPE == 3) then
+
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ b_displ_crust_mantle(i,1) = b_displ_crust_mantle(i,1) &
+ + b_deltat*b_veloc_crust_mantle(i,1) + b_deltatsqover2*b_accel_crust_mantle(i,1)
+ b_veloc_crust_mantle(i,1) = b_veloc_crust_mantle(i,1) &
+ + b_deltatover2*b_accel_crust_mantle(i,1)
+ b_accel_crust_mantle(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
b_displ_crust_mantle(:,i) = b_displ_crust_mantle(:,i) &
+ b_deltat*b_veloc_crust_mantle(:,i) + b_deltatsqover2*b_accel_crust_mantle(:,i)
@@ -16,6 +26,8 @@
+ b_deltatover2*b_accel_crust_mantle(:,i)
b_accel_crust_mantle(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
! outer core
do i=1,NGLOB_OUTER_CORE
b_displ_outer_core(i) = b_displ_outer_core(i) &
@@ -24,7 +36,17 @@
+ b_deltatover2*b_accel_outer_core(i)
b_accel_outer_core(i) = 0._CUSTOM_REAL
enddo
+
! inner core
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ b_displ_inner_core(i,1) = b_displ_inner_core(i,1) &
+ + b_deltat*b_veloc_inner_core(i,1) + b_deltatsqover2*b_accel_inner_core(i,1)
+ b_veloc_inner_core(i,1) = b_veloc_inner_core(i,1) &
+ + b_deltatover2*b_accel_inner_core(i,1)
+ b_accel_inner_core(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_INNER_CORE
b_displ_inner_core(:,i) = b_displ_inner_core(:,i) &
+ b_deltat*b_veloc_inner_core(:,i) + b_deltatsqover2*b_accel_inner_core(:,i)
@@ -32,6 +54,8 @@
+ b_deltatover2*b_accel_inner_core(:,i)
b_accel_inner_core(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
endif ! SIMULATION_TYPE == 3
! compute the maximum of the norm of the displacement
@@ -73,7 +97,7 @@
call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -88,14 +112,15 @@
b_buffer_send_faces,b_buffer_received_faces,npoin2D_max_all_CM_IC, &
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -112,7 +137,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
endif
@@ -179,7 +204,7 @@
! assemble all the contributions between slices using MPI
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
if (SIMULATION_TYPE == 3) then
@@ -206,7 +231,7 @@
call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -221,15 +246,16 @@
b_buffer_send_faces,b_buffer_received_faces,npoin2D_max_all_CM_IC, &
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -246,7 +272,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
do while (b_iphase <= 7) ! make sure the last communications are finished and processed
@@ -265,7 +291,7 @@
NGLOB2DMAX_XY_OC_VAL,NCHUNKS_VAL,b_iphase)
enddo
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
! Newmark time scheme - corrector for fluid parts
do i=1,NGLOB_OUTER_CORE
@@ -315,7 +341,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -326,11 +352,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -367,11 +393,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
endif
@@ -400,16 +426,16 @@
npoin2D_cube_from_slices,b_buffer_all_cube_from_slices,b_buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -438,19 +464,18 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
endif
! Stacey
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
if(SIMULATION_TYPE == 3) then
-
call compute_stacey_crust_mantle_backward(ichunk, &
NSTEP,it,ibool_crust_mantle, &
b_accel_crust_mantle, &
@@ -534,7 +559,7 @@
! crust/mantle and inner core handled in the same call
! in order to reduce the number of MPI messages by 2
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
if (SIMULATION_TYPE == 3) then
@@ -595,7 +620,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -606,11 +631,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -647,11 +672,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
! Deville routine
@@ -679,16 +704,16 @@
npoin2D_cube_from_slices,b_buffer_all_cube_from_slices,b_buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -717,12 +742,12 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
! assemble all the contributions between slices using MPI
@@ -759,10 +784,9 @@
enddo
endif ! end of assembling forces with the central cube
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
do i=1,NGLOB_CRUST_MANTLE
b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
@@ -770,9 +794,7 @@
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
enddo
-
else
-
do i=1,NGLOB_CRUST_MANTLE
b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
@@ -780,7 +802,6 @@
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
enddo
-
endif
endif ! SIMULATION_TYPE == 3
@@ -792,7 +813,7 @@
rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
ABSORBING_CONDITIONS)
endif
@@ -808,10 +829,18 @@
! Newmark time scheme - corrector for elastic parts
if (SIMULATION_TYPE == 3) then
+
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ b_veloc_crust_mantle(i,1) = b_veloc_crust_mantle(i,1) + b_deltatover2*b_accel_crust_mantle(i,1)
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
b_veloc_crust_mantle(:,i) = b_veloc_crust_mantle(:,i) + b_deltatover2*b_accel_crust_mantle(:,i)
enddo
+ endif
+
! inner core
do i=1,NGLOB_INNER_CORE
b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*rmass_inner_core(i) &
@@ -819,8 +848,18 @@
b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*rmass_inner_core(i) &
- b_two_omega_earth*b_veloc_inner_core(1,i)
b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*rmass_inner_core(i)
+ enddo
+
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ b_veloc_inner_core(i,1) = b_veloc_inner_core(i,1) + b_deltatover2*b_accel_inner_core(i,1)
+ enddo
+ else
+ do i=1,NGLOB_INNER_CORE
b_veloc_inner_core(:,i) = b_veloc_inner_core(:,i) + b_deltatover2*b_accel_inner_core(:,i)
enddo
+ endif
+
endif ! SIMULATION_TYPE == 3
! restores last time snapshot saved for backward/reconstruction of wavefields
@@ -874,13 +913,3 @@
seismo_current = 0
endif
-!
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!
-
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_undo_att.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part2_undo_att.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -8,7 +8,17 @@
! backward field
if (SIMULATION_TYPE == 3) then
+
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ b_displ_crust_mantle(i,1) = b_displ_crust_mantle(i,1) &
+ + b_deltat*b_veloc_crust_mantle(i,1) + b_deltatsqover2*b_accel_crust_mantle(i,1)
+ b_veloc_crust_mantle(i,1) = b_veloc_crust_mantle(i,1) &
+ + b_deltatover2*b_accel_crust_mantle(i,1)
+ b_accel_crust_mantle(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
b_displ_crust_mantle(:,i) = b_displ_crust_mantle(:,i) &
+ b_deltat*b_veloc_crust_mantle(:,i) + b_deltatsqover2*b_accel_crust_mantle(:,i)
@@ -16,6 +26,8 @@
+ b_deltatover2*b_accel_crust_mantle(:,i)
b_accel_crust_mantle(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
! outer core
do i=1,NGLOB_OUTER_CORE
b_displ_outer_core(i) = b_displ_outer_core(i) &
@@ -24,7 +36,17 @@
+ b_deltatover2*b_accel_outer_core(i)
b_accel_outer_core(i) = 0._CUSTOM_REAL
enddo
+
! inner core
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ b_displ_inner_core(i,1) = b_displ_inner_core(i,1) &
+ + b_deltat*b_veloc_inner_core(i,1) + b_deltatsqover2*b_accel_inner_core(i,1)
+ b_veloc_inner_core(i,1) = b_veloc_inner_core(i,1) &
+ + b_deltatover2*b_accel_inner_core(i,1)
+ b_accel_inner_core(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
do i=1,NGLOB_INNER_CORE
b_displ_inner_core(:,i) = b_displ_inner_core(:,i) &
+ b_deltat*b_veloc_inner_core(:,i) + b_deltatsqover2*b_accel_inner_core(:,i)
@@ -32,6 +54,8 @@
+ b_deltatover2*b_accel_inner_core(:,i)
b_accel_inner_core(:,i) = 0._CUSTOM_REAL
enddo
+ endif
+
endif ! SIMULATION_TYPE == 3
! compute the maximum of the norm of the displacement
@@ -77,7 +101,7 @@
call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -92,14 +116,15 @@
b_buffer_send_faces,b_buffer_received_faces,npoin2D_max_all_CM_IC, &
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -116,7 +141,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
endif
@@ -210,7 +235,7 @@
! assemble all the contributions between slices using MPI
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
if (SIMULATION_TYPE == 3) then
@@ -237,7 +262,7 @@
call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -252,15 +277,16 @@
b_buffer_send_faces,b_buffer_received_faces,npoin2D_max_all_CM_IC, &
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE, &
+ wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid, &
- b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core_dummy, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -277,7 +303,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core,MOVIE_VOLUME,&
- istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk,USE_LDDRK,SIMULATION_TYPE)
endif
do while (b_iphase <= 7) ! make sure the last communications are finished and processed
@@ -296,7 +322,7 @@
NGLOB2DMAX_XY_OC_VAL,NCHUNKS_VAL,b_iphase)
enddo
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
! Newmark time scheme - corrector for fluid parts
do i=1,NGLOB_OUTER_CORE
@@ -346,7 +372,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -357,11 +383,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -398,11 +424,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
endif
@@ -431,16 +457,16 @@
npoin2D_cube_from_slices,b_buffer_all_cube_from_slices,b_buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -469,12 +495,12 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
endif
@@ -599,7 +625,7 @@
! crust/mantle and inner core handled in the same call
! in order to reduce the number of MPI messages by 2
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
if (SIMULATION_TYPE == 3) then
@@ -660,7 +686,7 @@
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
muhstore_crust_mantle,eta_anisostore_crust_mantle, &
c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
@@ -671,11 +697,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -712,11 +738,11 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_crust_mantle,5),&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
endif
! Deville routine
@@ -744,16 +770,16 @@
npoin2D_cube_from_slices,b_buffer_all_cube_from_slices,b_buffer_slices,ibool_central_cube, &
receiver_cube_from_slices,ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,INCLUDE_CENTRAL_CUBE,b_iphase_CC, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core,wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -782,12 +808,12 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
- istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
+ size(factor_common_inner_core,5),&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
endif
! assemble all the contributions between slices using MPI
@@ -824,28 +850,44 @@
enddo
endif ! end of assembling forces with the central cube
-! ------------------- new non blocking implementation -------------------
+! ------------------- non blocking implementation -------------------
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
- do i=1,NGLOB_CRUST_MANTLE
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .and. .not. USE_LDDRK) then
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*b_rmassx_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
- b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*b_rmassy_crust_mantle(i) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
- b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
- enddo
-
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*b_rmassz_crust_mantle(i)
+ enddo
+ else
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i)
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i)
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
else
-
- do i=1,NGLOB_CRUST_MANTLE
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. .not. USE_LDDRK)then
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*b_rmassx_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
- b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*b_rmassy_crust_mantle(i) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
- b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
- enddo
-
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*b_rmassz_crust_mantle(i)
+ enddo
+ else
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*b_rmassz_crust_mantle(i) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i)
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*b_rmassz_crust_mantle(i) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i)
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*b_rmassz_crust_mantle(i)
+ enddo
+ endif
endif
endif ! SIMULATION_TYPE == 3
@@ -854,12 +896,12 @@
if(OCEANS_VAL) then
if(SIMULATION_TYPE == 3) then
call compute_coupling_ocean(b_accel_crust_mantle, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,b_rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
endif
endif
@@ -873,19 +915,47 @@
! Newmark time scheme - corrector for elastic parts
if (SIMULATION_TYPE == 3) then
+
! mantle
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_CRUST_MANTLE*NDIM
+ b_veloc_crust_mantle(i,1) = b_veloc_crust_mantle(i,1) + b_deltatover2*b_accel_crust_mantle(i,1)
+ enddo
+ else
do i=1,NGLOB_CRUST_MANTLE
b_veloc_crust_mantle(:,i) = b_veloc_crust_mantle(:,i) + b_deltatover2*b_accel_crust_mantle(:,i)
enddo
+ endif
+
! inner core
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. .not. USE_LDDRK) then
+ do i=1,NGLOB_INNER_CORE
+ b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*b_rmassx_inner_core(i) &
+ + b_two_omega_earth*b_veloc_inner_core(2,i)
+ b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*b_rmassy_inner_core(i) &
+ - b_two_omega_earth*b_veloc_inner_core(1,i)
+ b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*b_rmass_inner_core(i)
+ enddo
+ else
+ do i=1,NGLOB_INNER_CORE
+ b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*b_rmass_inner_core(i) &
+ + b_two_omega_earth*b_veloc_inner_core(2,i)
+ b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*b_rmass_inner_core(i) &
+ - b_two_omega_earth*b_veloc_inner_core(1,i)
+ b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*b_rmass_inner_core(i)
+ enddo
+ endif
+
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,NGLOB_INNER_CORE*NDIM
+ b_veloc_inner_core(i,1) = b_veloc_inner_core(i,1) + b_deltatover2*b_accel_inner_core(i,1)
+ enddo
+ else
do i=1,NGLOB_INNER_CORE
- b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*rmass_inner_core(i) &
- + b_two_omega_earth*b_veloc_inner_core(2,i)
- b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*rmass_inner_core(i) &
- - b_two_omega_earth*b_veloc_inner_core(1,i)
- b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*rmass_inner_core(i)
b_veloc_inner_core(:,i) = b_veloc_inner_core(:,i) + b_deltatover2*b_accel_inner_core(:,i)
enddo
+ endif
+
endif ! SIMULATION_TYPE == 3
! restores last time snapshot saved for backward/reconstruction of wavefields
@@ -914,7 +984,6 @@
ispec_selected_rec,number_receiver_global, &
seismo_current,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
seismograms)
-
endif
endif ! nrec_local
@@ -960,13 +1029,3 @@
seismo_current = 0
endif
-!
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!
-
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part3_kernel_computation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part3_kernel_computation.f90 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/part3_kernel_computation.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -8,25 +8,25 @@
rho_kl_crust_mantle,beta_kl_crust_mantle, &
alpha_kl_crust_mantle,cijkl_kl_crust_mantle, &
accel_crust_mantle,b_displ_crust_mantle, &
- deltat,displ_crust_mantle,hprime_xx,hprime_xxT,&
+ deltat,hprime_xx,hprime_xxT,&
xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle)
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,ANISOTROPIC_KL,&
+ epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
! outer core
call compute_kernels_outer_core(ibool_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
- hprime_xx,hprime_yy,hprime_zz, &
+ hprime_xx,hprime_xxT, &
displ_outer_core,accel_outer_core, &
b_displ_outer_core,b_accel_outer_core, &
vector_accel_outer_core,vector_displ_outer_core, &
b_vector_displ_outer_core, &
- div_displ_outer_core,b_div_displ_outer_core, &
+ div_displ_outer_core, &
rhostore_outer_core,kappavstore_outer_core, &
rho_kl_outer_core,alpha_kl_outer_core, &
- deviatoric_outercore,nspec_beta_kl_outer_core,beta_kl_outer_core, &
deltat)
! inner core
@@ -34,10 +34,11 @@
rho_kl_inner_core,beta_kl_inner_core, &
alpha_kl_inner_core, &
accel_inner_core,b_displ_inner_core, &
- deltat,displ_inner_core,hprime_xx,hprime_xxT,&
+ deltat,hprime_xx,hprime_xxT,&
xix_inner_core,xiy_inner_core,xiz_inner_core,&
etax_inner_core,etay_inner_core,etaz_inner_core,&
- gammax_inner_core,gammay_inner_core,gammaz_inner_core)
+ gammax_inner_core,gammay_inner_core,gammaz_inner_core,&
+ epsilondev_inner_core,eps_trace_over_3_inner_core)
! NOISE TOMOGRAPHY --- source strength kernel
if (NOISE_TOMOGRAPHY == 3) &
@@ -239,7 +240,7 @@
icb_kl = icb_kl + (icb_kl_top - icb_kl_bot) * deltat
endif
- ! approximate hessian
+ ! approximate Hessian
if( APPROXIMATE_HESS_KL ) then
call compute_kernels_hessian(ibool_crust_mantle, &
hess_kl_crust_mantle,&
@@ -249,10 +250,3 @@
endif ! end of if computing kernels
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------------------------
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk 2013-07-24 13:50:01 UTC (rev 22657)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk 2013-07-24 14:02:17 UTC (rev 22658)
@@ -98,7 +98,7 @@
# These files come from the shared directory
specfem3D_SHARED_OBJECTS = \
$O/auto_ner.o \
- $O/broadcast_compute_parameters.o \
+ $O/broadcast_computed_parameters.o \
$O/calendar.o \
$O/count_number_of_sources.o \
$O/create_name_database.o \
@@ -174,8 +174,8 @@
$O/compute_forces_outer_core_noDev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_outer_core_noDev.f90
${FCCOMPILE_CHECK} -c -o $O/compute_forces_outer_core_noDev.o ${FCFLAGS_f90} $S/compute_forces_outer_core_noDev.f90
-$O/compute_forces_outer_core_Dev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_outer_core_Dev.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_forces_outer_core_Dev.o ${FCFLAGS_f90} $S/compute_forces_outer_core_Dev.f90
+$O/compute_forces_outer_core_Dev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_outer_core_Dev.F90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_forces_outer_core_Dev.o ${FCFLAGS_f90} $S/compute_forces_outer_core_Dev.F90
$O/compute_forces_inner_core_noDev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_inner_core_noDev.f90
${FCCOMPILE_CHECK} -c -o $O/compute_forces_inner_core_noDev.o ${FCFLAGS_f90} $S/compute_forces_inner_core_noDev.f90
@@ -183,8 +183,8 @@
$O/compute_forces_inner_core_Dev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_inner_core_Dev.F90
${FCCOMPILE_CHECK} -c -o $O/compute_forces_inner_core_Dev.o ${FCFLAGS_f90} $S/compute_forces_inner_core_Dev.F90
-$O/compute_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_kernels.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o ${FCFLAGS_f90} $S/compute_kernels.f90
+$O/compute_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_kernels.F90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o ${FCFLAGS_f90} $S/compute_kernels.F90
$O/compute_seismograms.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_seismograms.f90
${FCCOMPILE_CHECK} -c -o $O/compute_seismograms.o ${FCFLAGS_f90} $S/compute_seismograms.f90
@@ -231,7 +231,7 @@
$O/define_derivation_matrices.o: ${SETUP}/constants.h $S/define_derivation_matrices.f90
${FCCOMPILE_CHECK} -c -o $O/define_derivation_matrices.o ${FCFLAGS_f90} $S/define_derivation_matrices.f90
-$O/get_attenuation.o: ${SETUP}/constants.h $S/get_attenuation.f90
+$O/get_attenuation.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/get_attenuation.f90
${FCCOMPILE_CHECK} -c -o $O/get_attenuation.o ${FCFLAGS_f90} $S/get_attenuation.f90
$O/get_backazimuth.o: ${SETUP}/constants.h $S/get_backazimuth.f90
@@ -315,7 +315,7 @@
$O/locate_receivers.o: ${SETUP}/constants.h $S/locate_receivers.f90
${MPIFCCOMPILE_CHECK} -c -o $O/locate_receivers.o ${FCFLAGS_f90} $S/locate_receivers.f90
-$O/locate_regular_points.o: ${SETUP}/constants.h $S/locate_regular_points.f90
+$O/locate_regular_points.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/locate_regular_points.f90
${MPIFCCOMPILE_CHECK} -c -o $O/locate_regular_points.o ${FCFLAGS_f90} $S/locate_regular_points.f90
$O/locate_sources.o: ${SETUP}/constants.h $S/locate_sources.f90
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/expand_compute_strain_product.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/expand_compute_strain_product.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/expand_compute_strain_product.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -0,0 +1,25 @@
+
+ program expand_compute_strain_product
+
+!! DK DK July 2013: expand the calculations in compute_strain_product() in order to speed up the code
+
+ implicit none
+
+ integer :: i,j,p
+
+ ! Computing the 21 strain products without assuming eps(i)*b_eps(j) = eps(j)*b_eps(i)
+ p=1
+ do i=1,6
+ do j=i,6
+ print *,'prod(',p,')=eps',i,'*b_eps',j
+ if(j>i) then
+ print *,'prod(',p,')=prod(',p,')+eps',j,'*b_eps',i
+ if(j>3 .and. i<4) print *,'prod(',p,') = prod(',p,') * 2._CUSTOM_REAL'
+ endif
+ if(i>3) print *,'prod(',p,') = prod(',p,') * 4._CUSTOM_REAL'
+ p=p+1
+ enddo
+ enddo
+
+ end program expand_compute_strain_product
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/test_shape_for_vectorization_if_not_first_index.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/test_shape_for_vectorization_if_not_first_index.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/test_shape_for_vectorization_if_not_first_index.f90 2013-07-24 14:02:17 UTC (rev 22658)
@@ -0,0 +1,64 @@
+
+program test_shape_vectorization
+
+!! DK DK July 2013: test if we can use (ijk,1,1) instead of (i,j,k) when the group is not located at the first index
+!! DK DK July 2013: (do NOT compile with range checking of course)
+
+implicit none
+
+integer, parameter :: NX = 8
+integer, parameter :: NY = 8
+integer, parameter :: NZ = 5
+integer, parameter :: NL = 7
+integer, parameter :: NM = 4
+
+integer, dimension(NL,NX,NY,NZ,NM) :: a
+
+integer :: i,j,k,l,m,ijk
+integer :: a_from_first_order,difference
+
+real :: r
+
+call random_seed()
+
+open(unit=10,file='first_result.dat',status='unknown')
+do m = 1,NM
+do k = 1,NZ
+do j = 1,NY
+do i = 1,NX
+do l = 1,NL
+ call random_number(r)
+ a(l,i,j,k,m) = nint(r * 20000.) ! create test values
+ print *,a(l,i,j,k,m)
+ write(10,*) a(l,i,j,k,m)
+enddo
+enddo
+enddo
+enddo
+enddo
+close(10)
+
+print *
+!! DK DK in practice it gives the exact same order, thus the trick works fine
+print *,'now in the vectorized order order'
+print *
+
+open(unit=10,file='first_result.dat',status='old')
+do m = 1,NM
+do ijk = 1,NX*NY*NZ
+do l = 1,NL
+ read(10,*) a_from_first_order
+ difference = abs(a(l,ijk,1,1,m) - a_from_first_order)
+ print *,a(l,ijk,1,1,m),difference
+ if(difference /= 0) stop 'error, difference between the two orders is not zero'
+enddo
+enddo
+enddo
+close(10)
+
+print *
+print *,'the test is successful, the two orders are 100% identical'
+print *
+
+end program test_shape_vectorization
+
More information about the CIG-COMMITS
mailing list