[cig-commits] r22735 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: . src/cuda src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Aug 29 00:47:36 PDT 2013
Author: danielpeter
Date: 2013-08-29 00:47:36 -0700 (Thu, 29 Aug 2013)
New Revision: 22735
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element_strain.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/broadcast_computed_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_attenuation_adios.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
Log:
adds new files compute_element_strain.F90 and update_displacement_Newmark.f90; prepares some more undoatt routines; updates restart file handling for restarting/checkpointing capabilities
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in 2013-08-29 07:47:36 UTC (rev 22735)
@@ -230,13 +230,13 @@
ifdef CLEAN
clean:
@echo "cleaning by CLEAN defined"
+ -rm -f $(foreach dir, $(CLEAN), $($(dir)_OBJECTS) $($(dir)_MODULES) $($(dir)_SHARED_OBJECTS) $($(dir)_TARGETS))
(rm -rf bin obj )
- -rm -f $(foreach dir, $(CLEAN), $($(dir)_OBJECTS) $($(dir)_MODULES) $($(dir)_SHARED_OBJECTS) $($(dir)_TARGETS))
else
clean:
@echo "cleaning by CLEAN not defined"
+ -rm -f $(foreach dir, $(SUBDIRS), $($(dir)_OBJECTS) $($(dir)_MODULES) $($(dir)_TARGETS))
(rm -rf bin obj )
- -rm -f $(foreach dir, $(SUBDIRS), $($(dir)_OBJECTS) $($(dir)_MODULES) $($(dir)_TARGETS))
endif
help:
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2013-08-29 07:47:36 UTC (rev 22735)
@@ -543,6 +543,7 @@
int absorbing_conditions;
int attenuation;
+ int undo_attenuation;
int use_attenuation_mimic;
int use_3d_attenuation_arrays;
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2013-08-29 07:47:36 UTC (rev 22735)
@@ -135,7 +135,7 @@
int* OCEANS_f,
int* GRAVITY_f,
int* ROTATION_f,
- int* ATTENUATION_f,
+ int* ATTENUATION_f,int* UNDO_ATTENUATION_f,
int* USE_ATTENUATION_MIMIC_f,
int* USE_3D_ATTENUATION_ARRAYS_f,
int* COMPUTE_AND_STORE_STRAIN_f,
@@ -238,6 +238,7 @@
mp->rotation = *ROTATION_f;
mp->attenuation = *ATTENUATION_f;
+ mp->undo_attenuation = *UNDO_ATTENUATION_f;
mp->use_attenuation_mimic = *USE_ATTENUATION_MIMIC_f;
mp->use_3d_attenuation_arrays = *USE_3D_ATTENUATION_ARRAYS_f;
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2013-08-29 07:47:36 UTC (rev 22735)
@@ -366,7 +366,7 @@
int* OCEANS_f,
int* GRAVITY_f,
int* ROTATION_f,
- int* ATTENUATION_f,
+ int* ATTENUATION_f,int* UNDO_ATTENUATION_f,
int* USE_ATTENUATION_MIMIC_f,
int* USE_3D_ATTENUATION_ARRAYS_f,
int* COMPUTE_AND_STORE_STRAIN_f,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/broadcast_computed_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/broadcast_computed_parameters.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/broadcast_computed_parameters.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -85,8 +85,10 @@
HONOR_1D_SPHERICAL_MOHO,MOVIE_COARSE, &
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY,&
ROTATE_SEISMOGRAMS_RT,WRITE_SEISMOGRAMS_BY_MASTER,USE_BINARY_FOR_LARGE_FILE, &
- SAVE_REGULAR_KL,PARTIAL_PHYS_DISPERSION_ONLY,UNDO_ATTENUATION, &
- USE_LDDRK,INCREASE_CFL_FOR_LDDRK,ANISOTROPIC_KL,SAVE_TRANSVERSE_KL_ONLY,APPROXIMATE_HESS_KL, &
+ SAVE_REGULAR_KL, &
+ PARTIAL_PHYS_DISPERSION_ONLY,UNDO_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, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -28,7 +28,14 @@
subroutine read_compute_parameters()
- use constants
+ use constants,only: &
+ TINYVAL,R_EARTH_KM,DEGREES_TO_RADIANS, &
+ SUPPRESS_CRUSTAL_MESH,ADD_4TH_DOUBLING, &
+ DO_BENCHMARK_RUN_ONLY,NSTEP_FOR_BENCHMARK, &
+ IREGION_CRUST_MANTLE,IREGION_INNER_CORE, &
+ NGLLX,NGLLY,NGLLZ
+
+
use shared_parameters
implicit none
@@ -113,14 +120,18 @@
! initial guess : compute total number of time steps, rounded to next multiple of 100
NSTEP = 100 * (int(RECORD_LENGTH_IN_MINUTES * 60.d0 / (100.d0*DT)) + 1)
+ ! if doing benchmark runs to measure scaling of the code for a limited number of time steps only
+ if (DO_BENCHMARK_RUN_ONLY) NSTEP = NSTEP_FOR_BENCHMARK
+
+ ! noise simulations
+ if ( NOISE_TOMOGRAPHY /= 0 ) then
+ ! time steps needs to be doubled, due to +/- branches
+ NSTEP = 2*NSTEP-1
+ endif
+
!! DK DK make sure NSTEP is a multiple of NT_DUMP_ATTENUATION
if(UNDO_ATTENUATION .and. mod(NSTEP,NT_DUMP_ATTENUATION) /= 0) NSTEP = (NSTEP/NT_DUMP_ATTENUATION + 1)*NT_DUMP_ATTENUATION
-! if doing benchmark runs to measure scaling of the code for a limited number of time steps only
- if (DO_BENCHMARK_RUN_ONLY) NSTEP = NSTEP_FOR_BENCHMARK
-
- if ( NOISE_TOMOGRAPHY /= 0 ) NSTEP = 2*NSTEP-1 ! time steps needs to be doubled, due to +/- branches
-
! subsets used to save seismograms must not be larger than the whole time series,
! otherwise we waste memory
if(NTSTEP_BETWEEN_OUTPUT_SEISMOS > NSTEP) NTSTEP_BETWEEN_OUTPUT_SEISMOS = NSTEP
@@ -132,11 +143,7 @@
240.d0/NEX_ETA*18.d0*ANGULAR_WIDTH_ETA_IN_DEGREES/90.d0)
! checks parameters
- call rcp_check_parameters(NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
- NCHUNKS,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
- ATTENUATION,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,OUTPUT_SEISMOS_SAC_ALPHANUM)
+ call rcp_check_parameters()
! check that mesh can be coarsened in depth three or four times
CUT_SUPERBRICK_XI=.false.
@@ -259,23 +266,16 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine rcp_check_parameters(NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
- NCHUNKS,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
- ATTENUATION,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,OUTPUT_SEISMOS_SAC_ALPHANUM)
+ subroutine rcp_check_parameters()
- use constants
+ use constants,only: &
+ CUSTOM_REAL,SIZE_REAL,SIZE_DOUBLE,NUMFACES_SHARED,NUMCORNERS_SHARED, &
+ N_SLS,NGNOD,NGNOD2D
+ use shared_parameters
+
implicit none
- integer NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NCHUNKS,NTSTEP_BETWEEN_OUTPUT_SEISMOS
-
- double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES
-
- logical ATTENUATION,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,OUTPUT_SEISMOS_SAC_ALPHANUM
-
! checks parameters
if(NCHUNKS /= 1 .and. NCHUNKS /= 2 .and. NCHUNKS /= 3 .and. NCHUNKS /= 6) &
@@ -295,6 +295,21 @@
if(ABSORBING_CONDITIONS .and. NCHUNKS == 3) &
stop 'absorbing conditions not supported for three chunks yet'
+ if(ATTENUATION_3D .and. .not. ATTENUATION) &
+ stop 'need ATTENUATION to use ATTENUATION_3D'
+
+ if(SAVE_TRANSVERSE_KL_ONLY .and. .not. ANISOTROPIC_KL) &
+ stop 'need ANISOTROPIC_KL to use SAVE_TRANSVERSE_KL_ONLY'
+
+ if(PARTIAL_PHYS_DISPERSION_ONLY .and. UNDO_ATTENUATION) &
+ stop 'cannot have both PARTIAL_PHYS_DISPERSION_ONLY and UNDO_ATTENUATION, they are mutually exclusive'
+
+ if(UNDO_ATTENUATION .and. MOVIE_VOLUME .and. MOVIE_VOLUME_TYPE == 4 ) &
+ stop 'MOVIE_VOLUME_TYPE == 4 is not implemented for UNDO_ATTENUATION in order to save memory'
+
+ if(UNDO_ATTENUATION .and. NUMBER_OF_THIS_RUN > 1) &
+ stop 'we currently do not support NUMBER_OF_THIS_RUN > 1 in the case of UNDO_ATTENUATION'
+
if (OUTPUT_SEISMOS_SAC_ALPHANUM .and. (mod(NTSTEP_BETWEEN_OUTPUT_SEISMOS,5)/=0)) &
stop 'if OUTPUT_SEISMOS_SAC_ALPHANUM = .true. then NTSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of 5, check the Par_file'
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -247,11 +247,23 @@
stop 'an error occurred while reading the parameter file'
endif
+!daniel debug: status of implementation
+
!! DK DK July 2013: temporary, the time for Matthieu Lefebvre to merge his ADIOS implementation
if( ADIOS_ENABLED ) then
- stop 'ADIOS support not implemented yet'
+ stop 'ADIOS_ENABLED support not implemented yet'
endif
+ if( USE_LDDRK ) then
+ stop 'USE_LDDRK support not implemented yet'
+ endif
+ if( EXACT_MASS_MATRIX_FOR_ROTATION ) then
+ stop 'EXACT_MASS_MATRIX_FOR_ROTATION support not implemented yet'
+ endif
+ if( UNDO_ATTENUATION ) then
+ stop 'UNDO_ATTENUATION support not implemented yet'
+ endif
+
end subroutine read_parameter_file
!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -38,8 +38,8 @@
use specfem_par,only: &
GPU_MODE,Mesh_pointer, &
COMPUTE_AND_STORE_STRAIN, &
- SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0, &
- NSTEP,it, &
+ SIMULATION_TYPE,time_start,DT,t0, &
+ NSTEP,it,it_begin,it_end,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN, &
myrank
use specfem_par_crustmantle,only: displ_crust_mantle,b_displ_crust_mantle, &
@@ -63,11 +63,14 @@
integer :: ihours,iminutes,iseconds,int_tCPU, &
ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
ihours_total,iminutes_total,iseconds_total,int_t_total
+ integer :: it_run,nstep_run
+ logical :: SHOW_SEPARATE_RUN_INFORMATION
! to determine date and time at which the run will finish
character(len=8) datein
character(len=10) timein
character(len=5) :: zone
integer, dimension(8) :: time_values
+
character(len=3), dimension(12) :: month_name
character(len=3), dimension(0:6) :: weekday_name
data month_name /'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
@@ -159,6 +162,38 @@
if(myrank == 0) then
+ ! this is in the case of restart files, when a given run consists of several partial runs
+ ! information about the current run only
+ SHOW_SEPARATE_RUN_INFORMATION = ( NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS )
+ it_run = it - it_begin + 1
+ nstep_run = it_end - it_begin + 1
+
+ ! elapsed time since beginning of the simulation
+ tCPU = wtime() - time_start
+ int_tCPU = int(tCPU)
+ ihours = int_tCPU / 3600
+ iminutes = (int_tCPU - 3600*ihours) / 60
+ iseconds = int_tCPU - 3600*ihours - 60*iminutes
+
+ ! compute estimated remaining simulation time
+ if (SHOW_SEPARATE_RUN_INFORMATION) then
+ t_remain = (it_end - it) * (tCPU/dble(it_run))
+ else
+ t_remain = (NSTEP - it) * (tCPU/dble(it))
+ endif
+ int_t_remain = int(t_remain)
+ ihours_remain = int_t_remain / 3600
+ iminutes_remain = (int_t_remain - 3600*ihours_remain) / 60
+ iseconds_remain = int_t_remain - 3600*ihours_remain - 60*iminutes_remain
+
+ ! compute estimated total simulation time
+ t_total = t_remain + tCPU
+ int_t_total = int(t_total)
+ ihours_total = int_t_total / 3600
+ iminutes_total = (int_t_total - 3600*ihours_total) / 60
+ iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
+
+ ! user output
write(IMAIN,*) 'Time step # ',it
write(IMAIN,*) 'Time: ',sngl(((it-1)*DT-t0)/60.d0),' minutes'
@@ -178,39 +213,31 @@
write(IMAIN,*) 'Max of strain, epsilondev_crust_mantle =',Strain2_norm_all
endif
- ! elapsed time since beginning of the simulation
- tCPU = wtime() - time_start
- int_tCPU = int(tCPU)
- ihours = int_tCPU / 3600
- iminutes = (int_tCPU - 3600*ihours) / 60
- iseconds = int_tCPU - 3600*ihours - 60*iminutes
write(IMAIN,*) 'Elapsed time in seconds = ',tCPU
write(IMAIN,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
- ! compute estimated remaining simulation time
- t_remain = (NSTEP - it) * (tCPU/dble(it))
- int_t_remain = int(t_remain)
- ihours_remain = int_t_remain / 3600
- iminutes_remain = (int_t_remain - 3600*ihours_remain) / 60
- iseconds_remain = int_t_remain - 3600*ihours_remain - 60*iminutes_remain
- write(IMAIN,*) 'Time steps done = ',it,' out of ',NSTEP
- write(IMAIN,*) 'Time steps remaining = ',NSTEP - it
+ if (SHOW_SEPARATE_RUN_INFORMATION) then
+ write(IMAIN,*) 'Time steps done for this run = ',it_run,' out of ',nstep_run
+ write(IMAIN,*) 'Time steps done in total = ',it,' out of ',NSTEP
+ write(IMAIN,*) 'Time steps remaining for this run = ',it_end - it
+ write(IMAIN,*) 'Time steps remaining for all runs = ',NSTEP - it
+ else
+ write(IMAIN,*) 'Time steps done = ',it,' out of ',NSTEP
+ write(IMAIN,*) 'Time steps remaining = ',NSTEP - it
+ endif
+
write(IMAIN,*) 'Estimated remaining time in seconds = ',t_remain
+ write(IMAIN,*) 'Estimated remaining time in seconds = ',t_remain
write(IMAIN,"(' Estimated remaining time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
ihours_remain,iminutes_remain,iseconds_remain
- ! compute estimated total simulation time
- t_total = t_remain + tCPU
- int_t_total = int(t_total)
- ihours_total = int_t_total / 3600
- iminutes_total = (int_t_total - 3600*ihours_total) / 60
- iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
write(IMAIN,*) 'Estimated total run time in seconds = ',t_total
write(IMAIN,"(' Estimated total run time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
ihours_total,iminutes_total,iseconds_total
write(IMAIN,*) 'We have done ',sngl(100.d0*dble(it)/dble(NSTEP)),'% of that'
+
if(it < NSTEP) then
! get current date
@@ -287,14 +314,12 @@
call flush_IMAIN()
! write time stamp file to give information about progression of simulation
- call write_timestamp_file(it,NSTEP,DT,t0,SIMULATION_TYPE, &
- Usolidnorm_all,Ufluidnorm_all,b_Usolidnorm_all,b_Ufluidnorm_all, &
+ call write_timestamp_file(Usolidnorm_all,Ufluidnorm_all,b_Usolidnorm_all,b_Ufluidnorm_all, &
tCPU,ihours,iminutes,iseconds, &
t_remain,ihours_remain,iminutes_remain,iseconds_remain, &
t_total,ihours_total,iminutes_total,iseconds_total, &
day_of_week,mon,day,year,hr,minutes, &
- day_of_week_remote,mon_remote,day_remote,year_remote,hr_remote,minutes_remote, &
- OUTPUT_FILES)
+ day_of_week_remote,mon_remote,day_remote,year_remote,hr_remote,minutes_remote)
! check stability of the code, exit if unstable
! negative values can occur with some compilers when the unstable value is greater
@@ -338,25 +363,22 @@
!------------------------------------------------------------------------------------------------------------------
!
- subroutine write_timestamp_file(it,NSTEP,DT,t0,SIMULATION_TYPE, &
- Usolidnorm_all,Ufluidnorm_all,b_Usolidnorm_all,b_Ufluidnorm_all, &
- tCPU,ihours,iminutes,iseconds, &
- t_remain,ihours_remain,iminutes_remain,iseconds_remain, &
- t_total,ihours_total,iminutes_total,iseconds_total, &
- day_of_week,mon,day,year,hr,minutes, &
- day_of_week_remote,mon_remote,day_remote,year_remote,hr_remote,minutes_remote, &
- OUTPUT_FILES)
+ subroutine write_timestamp_file(Usolidnorm_all,Ufluidnorm_all,b_Usolidnorm_all,b_Ufluidnorm_all, &
+ tCPU,ihours,iminutes,iseconds, &
+ t_remain,ihours_remain,iminutes_remain,iseconds_remain, &
+ t_total,ihours_total,iminutes_total,iseconds_total, &
+ day_of_week,mon,day,year,hr,minutes, &
+ day_of_week_remote,mon_remote,day_remote,year_remote,hr_remote,minutes_remote)
- use constants_solver
+ use constants,only: CUSTOM_REAL,IOUT, &
+ ADD_TIME_ESTIMATE_ELSEWHERE,HOURS_TIME_DIFFERENCE,MINUTES_TIME_DIFFERENCE
+ use specfem_par,only: &
+ SIMULATION_TYPE,OUTPUT_FILES,DT,t0, &
+ NSTEP,it,it_begin,it_end,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN
+
implicit none
- ! time step
- integer :: it,NSTEP
- integer :: SIMULATION_TYPE
-
- double precision :: DT,t0
-
! maximum of the norm of the displacement and of the potential in the fluid
real(kind=CUSTOM_REAL) Usolidnorm_all,Ufluidnorm_all
real(kind=CUSTOM_REAL) b_Usolidnorm_all,b_Ufluidnorm_all
@@ -370,12 +392,11 @@
integer :: year,mon,day,hr,minutes,day_of_week, &
year_remote,mon_remote,day_remote,hr_remote,minutes_remote,day_of_week_remote
- character(len=150) :: OUTPUT_FILES
! local parameters
+ integer :: it_run,nstep_run
! names of the data files for all the processors in MPI
character(len=150) :: outputname
-
! to determine date and time at which the run will finish
character(len=3), dimension(12) :: month_name
character(len=3), dimension(0:6) :: weekday_name
@@ -405,8 +426,20 @@
write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
write(IOUT,*)
- write(IOUT,*) 'Time steps done = ',it,' out of ',NSTEP
- write(IOUT,*) 'Time steps remaining = ',NSTEP - it
+ if( NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS ) then
+ ! this is in the case of restart files, when a given run consists of several partial runs
+ ! information about the current run only
+ it_run = it - it_begin + 1
+ nstep_run = it_end - it_begin + 1
+ write(IOUT,*) 'Time steps done for this run = ',it_run,' out of ',nstep_run
+ write(IOUT,*) 'Time steps done in total = ',it,' out of ',NSTEP
+ write(IOUT,*) 'Time steps remaining for this run = ',it_end - it
+ write(IOUT,*) 'Time steps remaining for all runs = ',NSTEP - it
+ else
+ write(IOUT,*) 'Time steps done = ',it,' out of ',NSTEP
+ write(IOUT,*) 'Time steps remaining = ',NSTEP - it
+ endif
+
write(IOUT,*) 'Estimated remaining time in seconds = ',t_remain
write(IOUT,"(' Estimated remaining time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
ihours_remain,iminutes_remain,iseconds_remain
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -26,18 +26,18 @@
!=====================================================================
subroutine compute_element_iso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- kappavstore,muvstore, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc, &
- epsilondev_loc,rho_s_H,is_backward_field)
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ kappavstore,muvstore, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ epsilondev_loc,rho_s_H,is_backward_field)
use constants_solver
use specfem_par,only: COMPUTE_AND_STORE_STRAIN
@@ -64,7 +64,7 @@
kappavstore,muvstore
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
! attenuation
! memory variables for attenuation
@@ -74,7 +74,7 @@
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
- real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: one_minus_sum_beta
! gravity
double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
@@ -116,6 +116,12 @@
integer :: int_radius
integer :: iglob
+ logical :: dummyl
+
+!daniel debug
+ ! dummy to avoid compiler warning
+ dummyl = is_backward_field
+
! isotropic element
do k=1,NGLLZ
@@ -183,13 +189,14 @@
mul = muvstore(i,j,k,ispec)
! use unrelaxed parameters if attenuation
- if(ATTENUATION_VAL .and. (ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL)) then
+ if(ATTENUATION_VAL ) then
! precompute terms for attenuation if needed
- one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
mul = mul * one_minus_sum_beta_use
- else if( ATTENUATION_VAL ) then
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
- mul = mul * one_minus_sum_beta_use
endif
lambdalplus2mul = kappal + FOUR_THIRDS * mul
@@ -207,7 +214,7 @@
! subtract memory variables if attenuation
if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL ) then
-!daniel: att - debug update
+!daniel debug: att - debug update
! call compute_element_att_mem_up_cm(ispec,i,j,k, &
! R_xx(1,i,j,k,ispec), &
! R_yy(1,i,j,k,ispec), &
@@ -215,10 +222,9 @@
! R_xz(1,i,j,k,ispec), &
! R_yz(1,i,j,k,ispec), &
! epsilondev_loc(:,i,j,k),muvstore(i,j,k,ispec),is_backward_field)
-! dummy to avoid compiler warning
- if( is_backward_field ) then
- endif
+ ! note: function inlining is generally done by fortran compilers;
+ ! compilers decide based on performance heuristics
! note: fortran passes pointers to array location, thus R_memory(1,1,...) should be fine
call compute_element_att_stress(R_xx(1,i,j,k,ispec), &
R_yy(1,i,j,k,ispec), &
@@ -363,18 +369,18 @@
!
subroutine compute_element_tiso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc, &
- epsilondev_loc,rho_s_H,is_backward_field)
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ epsilondev_loc,rho_s_H,is_backward_field)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -413,10 +419,10 @@
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
- real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: one_minus_sum_beta
! gravity
double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
@@ -471,6 +477,12 @@
integer :: int_radius
integer :: iglob
+ logical :: dummyl
+
+!daniel debug
+ ! dummy to avoid compiler warning
+ dummyl = is_backward_field
+
! transverse isotropic element
do k=1,NGLLZ
@@ -546,15 +558,15 @@
! use unrelaxed parameters if attenuation
! eta does not need to be shifted since it is a ratio
- if(ATTENUATION_VAL .and. (ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL)) then
+ if( ATTENUATION_VAL ) then
! precompute terms for attenuation if needed
- one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
muvl = muvl * one_minus_sum_beta_use
muhl = muhl * one_minus_sum_beta_use
- else if( ATTENUATION_VAL ) then
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
- muvl = muvl * one_minus_sum_beta_use
- muhl = muhl * one_minus_sum_beta_use
endif
rhovpvsq = kappavl + FOUR_THIRDS * muvl !!! that is C
@@ -748,7 +760,9 @@
! subtract memory variables if attenuation
if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL ) then
- ! note: Fortran passes pointers to array location, thus R_memory(1,1,...) is fine
+ ! note: function inlining is generally done by fortran compilers;
+ ! compilers decide based on performance heuristics
+ ! note: fortran passes pointers to array location, thus R_memory(1,1,...) is fine
call compute_element_att_stress(R_xx(1,i,j,k,ispec), &
R_yy(1,i,j,k,ispec), &
R_xy(1,i,j,k,ispec), &
@@ -891,20 +905,20 @@
!
subroutine compute_element_aniso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
- c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
- c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc, &
- epsilondev_loc,rho_s_H,is_backward_field)
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ epsilondev_loc,rho_s_H,is_backward_field)
use constants_solver
use specfem_par,only: COMPUTE_AND_STORE_STRAIN
@@ -941,10 +955,10 @@
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
- real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: one_minus_sum_beta
! gravity
double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
@@ -988,6 +1002,12 @@
integer :: int_radius
integer :: iglob
+ logical :: dummyl
+
+!daniel debug
+ ! dummy to avoid compiler warning
+ dummyl = is_backward_field
+
! anisotropic elements
do k=1,NGLLZ
@@ -1072,10 +1092,13 @@
c56 = c56store(i,j,k,ispec)
c66 = c66store(i,j,k,ispec)
- if(ATTENUATION_VAL .and. (ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL)) then
+ if( ATTENUATION_VAL ) then
! precompute terms for attenuation if needed
- minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0_CUSTOM_REAL
- !mul = c44
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0_CUSTOM_REAL
+ else
+ minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0_CUSTOM_REAL
+ endif
mul = c44 * minus_sum_beta
c11 = c11 + FOUR_THIRDS * mul ! * minus_sum_beta * mul
c12 = c12 - TWO_THIRDS * mul
@@ -1086,19 +1109,6 @@
c44 = c44 + mul
c55 = c55 + mul
c66 = c66 + mul
- else if( ATTENUATION_VAL ) then
- minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0_CUSTOM_REAL
- !mul = c44
- mul = c44 * minus_sum_beta
- c11 = c11 + FOUR_THIRDS * mul ! * minus_sum_beta * mul
- c12 = c12 - TWO_THIRDS * mul
- c13 = c13 - TWO_THIRDS * mul
- c22 = c22 + FOUR_THIRDS * mul
- c23 = c23 - TWO_THIRDS * mul
- c33 = c33 + FOUR_THIRDS * mul
- c44 = c44 + mul
- c55 = c55 + mul
- c66 = c66 + mul
endif
sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
@@ -1121,6 +1131,8 @@
! subtract memory variables if attenuation
if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL ) then
+ ! note: function inlining is generally done by fortran compilers;
+ ! compilers decide based on performance heuristics
! note: Fortran passes pointers to array location, thus R_memory(1,1,...) is fine
call compute_element_att_stress(R_xx(1,i,j,k,ispec), &
@@ -1268,7 +1280,7 @@
subroutine compute_element_att_stress(R_xx_loc,R_yy_loc,R_xy_loc,R_xz_loc,R_yz_loc, &
sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz)
- use constants_solver
+ use constants_solver,only: CUSTOM_REAL,N_SLS
implicit none
@@ -1307,11 +1319,11 @@
!
subroutine compute_element_att_memory_cm(ispec,R_xx,R_yy,R_xy,R_xz,R_yz, &
- vx,vy,vz,vnspec,factor_common, &
- alphaval,betaval,gammaval, &
- c44store,muvstore, &
- epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
- epsilondev_loc,is_backward_field)
+ vnspec,factor_common, &
+ alphaval,betaval,gammaval, &
+ c44store,muvstore, &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
+ epsilondev_loc,is_backward_field)
! crust mantle
! update memory variables based upon the Runge-Kutta scheme
@@ -1343,9 +1355,9 @@
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUATION) :: R_xx,R_yy,R_xy,R_xz,R_yz
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
- real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: c44store
@@ -1409,11 +1421,11 @@
!
subroutine compute_element_att_memory_ic(ispec,R_xx,R_yy,R_xy,R_xz,R_yz, &
- vx,vy,vz,vnspec,factor_common, &
- alphaval,betaval,gammaval, &
- muvstore, &
- epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
- epsilondev_loc,is_backward_field)
+ vnspec,factor_common, &
+ alphaval,betaval,gammaval, &
+ muvstore, &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, &
+ epsilondev_loc,is_backward_field)
! inner core
! update memory variables based upon the Runge-Kutta scheme
@@ -1442,9 +1454,9 @@
R_xx,R_yy,R_xy,R_xz,R_yz
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
- real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: muvstore
@@ -1596,3 +1608,4 @@
enddo ! i_SLS
end subroutine compute_element_att_mem_up_cm
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element_strain.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element_strain.F90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element_strain.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -0,0 +1,675 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
+! --------------------------------------------------
+!
+! 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
+! August 2013
+!
+! 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.
+!
+!=====================================================================
+
+
+ subroutine compute_element_strain_undo_att_Dev(ispec,nglob,nspec, &
+ displ,ibool, &
+ hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ epsilondev_loc,eps_trace_over_3_loc)
+
+! computes strain for single element based on Deville routine setup (NGLLX == NGLLY == NGLLZ == 5)
+
+ use constants
+
+ implicit none
+
+ integer,intent(in) :: ispec,nglob,nspec
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob),intent(in) :: displ
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX),intent(in) :: hprime_xx,hprime_xxT
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ),intent(out) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ),intent(out) :: eps_trace_over_3_loc
+
+! local variable
+ integer :: i,j,k,iglob
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(dummyy_loc,B2_m1_m2_5points)
+ equivalence(dummyz_loc,B3_m1_m2_5points)
+
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(tempy1,C2_m1_m2_5points)
+ equivalence(tempz1,C3_m1_m2_5points)
+
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: &
+ C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(tempy3,C2_mxm_m2_m1_5points)
+ equivalence(tempz3,C3_mxm_m2_m1_5points)
+
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: &
+ A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+ equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duydyl,duzdzl,duxdyl,duydxl,duzdxl,duxdzl,duzdyl,duydzl,&
+ duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl,&
+ duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#endif
+
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ iglob = ibool(ijk,1,1,ispec)
+ dummyx_loc(ijk,1,1) = displ(1,iglob)
+ dummyy_loc(ijk,1,1) = displ(2,iglob)
+ dummyz_loc(ijk,1,1) = displ(3,iglob)
+ enddo
+#else
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob)
+ enddo
+ enddo
+ enddo
+#endif
+
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+
+ C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+
+ C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(ijk,1,1,ispec)
+ xiyl = xiy(ijk,1,1,ispec)
+ xizl = xiz(ijk,1,1,ispec)
+ etaxl = etax(ijk,1,1,ispec)
+ etayl = etay(ijk,1,1,ispec)
+ etazl = etaz(ijk,1,1,ispec)
+ gammaxl = gammax(ijk,1,1,ispec)
+ gammayl = gammay(ijk,1,1,ispec)
+ gammazl = gammaz(ijk,1,1,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(ijk,1,1) + etaxl*tempx2(ijk,1,1) + gammaxl*tempx3(ijk,1,1)
+ duxdyl = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
+ duxdzl = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
+
+ duydxl = xixl*tempy1(ijk,1,1) + etaxl*tempy2(ijk,1,1) + gammaxl*tempy3(ijk,1,1)
+ duydyl = xiyl*tempy1(ijk,1,1) + etayl*tempy2(ijk,1,1) + gammayl*tempy3(ijk,1,1)
+ duydzl = xizl*tempy1(ijk,1,1) + etazl*tempy2(ijk,1,1) + gammazl*tempy3(ijk,1,1)
+
+ duzdxl = xixl*tempz1(ijk,1,1) + etaxl*tempz2(ijk,1,1) + gammaxl*tempz3(ijk,1,1)
+ duzdyl = xiyl*tempz1(ijk,1,1) + etayl*tempz2(ijk,1,1) + gammayl*tempz3(ijk,1,1)
+ duzdzl = xizl*tempz1(ijk,1,1) + etazl*tempz2(ijk,1,1) + gammazl*tempz3(ijk,1,1)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ ! strains
+ eps_trace_over_3_loc(ijk,1,1) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,ijk,1,1) = duxdxl - eps_trace_over_3_loc(ijk,1,1)
+ epsilondev_loc(2,ijk,1,1) = duydyl - eps_trace_over_3_loc(ijk,1,1)
+ epsilondev_loc(3,ijk,1,1) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc(4,ijk,1,1) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc(5,ijk,1,1) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ enddo
+#else
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ ! strains
+ eps_trace_over_3_loc(i,j,k) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,i,j,k) = duxdxl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(2,i,j,k) = duydyl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(3,i,j,k) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ enddo
+ enddo
+ enddo
+#endif
+
+ end subroutine compute_element_strain_undo_att_Dev
+
+
+!
+!--------------------------------------------------------------------------------------------
+!
+
+ subroutine compute_element_strain_undo_att_noDev(ispec,nglob,nspec, &
+ displ, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ ibool,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ epsilondev_loc,eps_trace_over_3_loc)
+
+! computes strain for single element
+
+ use constants
+
+ implicit none
+
+ ! element id
+ integer,intent(in) :: ispec
+
+ integer,intent(in) :: NSPEC,NGLOB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB),intent(in) :: displ
+
+ ! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX),intent(in) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY),intent(in) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ),intent(in) :: hprime_zz
+
+ ! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC),intent(in) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC),intent(in) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ),intent(out) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ),intent(out) :: eps_trace_over_3_loc
+
+ ! local parameters
+ integer :: iglob
+ integer :: i,j,k,l
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute deviatoric strain
+ eps_trace_over_3_loc(i,j,k) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,i,j,k) = duxdxl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(2,i,j,k) = duydyl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(3,i,j,k) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+
+ enddo ! NGLLX
+ enddo ! NGLLY
+ enddo ! NGLLZ
+
+ end subroutine compute_element_strain_undo_att_noDev
+
+
+!
+!--------------------------------------------------------------------------------------------
+!
+
+ subroutine compute_element_strain_att_Dev(ispec,nglob,nspec, &
+ displ,veloc,deltat,ibool, &
+ hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ epsilondev_loc_nplus1,eps_trace_over_3_loc_nplus1)
+
+ use constants
+
+ implicit none
+
+ integer :: ispec,nglob,nspec
+ real(kind=CUSTOM_REAL) :: deltat
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: displ,veloc
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_nplus1
+
+ ! local variable
+ integer :: i,j,k,iglob
+ real(kind=CUSTOM_REAL) :: templ
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(dummyy_loc,B2_m1_m2_5points)
+ equivalence(dummyz_loc,B3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(tempy1,C2_m1_m2_5points)
+ equivalence(tempz1,C3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: &
+ C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(tempy3,C2_mxm_m2_m1_5points)
+ equivalence(tempz3,C3_mxm_m2_m1_5points)
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: &
+ A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+ equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duydyl,duzdzl,duxdyl,duydxl,duzdxl,duxdzl,duzdyl,duydzl,&
+ duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl,&
+ duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#endif
+
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ iglob = ibool(ijk,1,1,ispec)
+ dummyx_loc(ijk,1,1) = displ(1,iglob) + deltat * veloc(1,iglob)
+ dummyy_loc(ijk,1,1) = displ(2,iglob) + deltat * veloc(2,iglob)
+ dummyz_loc(ijk,1,1) = displ(3,iglob) + deltat * veloc(3,iglob)
+ enddo
+#else
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob) + deltat * veloc(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob) + deltat * veloc(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob) + deltat * veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+#endif
+
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+
+ C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+
+ C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(ijk,1,1,ispec)
+ xiyl = xiy(ijk,1,1,ispec)
+ xizl = xiz(ijk,1,1,ispec)
+ etaxl = etax(ijk,1,1,ispec)
+ etayl = etay(ijk,1,1,ispec)
+ etazl = etaz(ijk,1,1,ispec)
+ gammaxl = gammax(ijk,1,1,ispec)
+ gammayl = gammay(ijk,1,1,ispec)
+ gammazl = gammaz(ijk,1,1,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(ijk,1,1) + etaxl*tempx2(ijk,1,1) + gammaxl*tempx3(ijk,1,1)
+ duxdyl = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
+ duxdzl = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
+
+ duydxl = xixl*tempy1(ijk,1,1) + etaxl*tempy2(ijk,1,1) + gammaxl*tempy3(ijk,1,1)
+ duydyl = xiyl*tempy1(ijk,1,1) + etayl*tempy2(ijk,1,1) + gammayl*tempy3(ijk,1,1)
+ duydzl = xizl*tempy1(ijk,1,1) + etazl*tempy2(ijk,1,1) + gammazl*tempy3(ijk,1,1)
+
+ duzdxl = xixl*tempz1(ijk,1,1) + etaxl*tempz2(ijk,1,1) + gammaxl*tempz3(ijk,1,1)
+ duzdyl = xiyl*tempz1(ijk,1,1) + etayl*tempz2(ijk,1,1) + gammayl*tempz3(ijk,1,1)
+ duzdzl = xizl*tempz1(ijk,1,1) + etazl*tempz2(ijk,1,1) + gammazl*tempz3(ijk,1,1)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ eps_trace_over_3_loc_nplus1 = templ
+ epsilondev_loc_nplus1(1,ijk,1,1) = duxdxl - templ
+ epsilondev_loc_nplus1(2,ijk,1,1) = duydyl - templ
+ epsilondev_loc_nplus1(3,ijk,1,1) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc_nplus1(4,ijk,1,1) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc_nplus1(5,ijk,1,1) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ enddo
+#else
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ eps_trace_over_3_loc_nplus1 = templ
+ epsilondev_loc_nplus1(1,i,j,k) = duxdxl - templ
+ epsilondev_loc_nplus1(2,i,j,k) = duydyl - templ
+ epsilondev_loc_nplus1(3,i,j,k) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc_nplus1(4,i,j,k) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc_nplus1(5,i,j,k) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ enddo
+ enddo
+ enddo
+#endif
+
+ end subroutine compute_element_strain_att_Dev
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -271,13 +271,13 @@
! (multiply by the inverse of the mass matrix and update velocity)
if(.NOT. GPU_MODE) then
! on CPU
- call compute_forces_ac_update_veloc(NGLOB_OUTER_CORE,veloc_outer_core,accel_outer_core, &
- deltatover2,rmass_outer_core)
+ call update_veloc_acoustic(NGLOB_OUTER_CORE,veloc_outer_core,accel_outer_core, &
+ deltatover2,rmass_outer_core)
! adjoint / kernel runs
if (SIMULATION_TYPE == 3) &
- call compute_forces_ac_update_veloc(NGLOB_OUTER_CORE_ADJOINT,b_veloc_outer_core,b_accel_outer_core, &
- b_deltatover2,rmass_outer_core)
+ call update_veloc_acoustic(NGLOB_OUTER_CORE_ADJOINT,b_veloc_outer_core,b_accel_outer_core, &
+ b_deltatover2,rmass_outer_core)
else
! on GPU
@@ -287,35 +287,4 @@
end subroutine compute_forces_acoustic
-!=====================================================================
- subroutine compute_forces_ac_update_veloc(NGLOB,veloc_outer_core,accel_outer_core, &
- deltatover2,rmass_outer_core)
-
- use constants_solver,only: CUSTOM_REAL
-
- implicit none
-
- integer :: NGLOB
-
- ! velocity potential
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: veloc_outer_core,accel_outer_core
-
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmass_outer_core
-
- real(kind=CUSTOM_REAL) :: deltatover2
-
- ! local parameters
- integer :: i
-
- ! Newmark time scheme
- ! multiply by the inverse of the mass matrix and update velocity
-
- do i=1,NGLOB
- accel_outer_core(i) = accel_outer_core(i)*rmass_outer_core(i)
- veloc_outer_core(i) = veloc_outer_core(i) + deltatover2*accel_outer_core(i)
- enddo
-
- end subroutine compute_forces_ac_update_veloc
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -36,7 +36,7 @@
epsilondev_xz,epsilondev_yz, &
epsilon_trace_over_3, &
alphaval,betaval,gammaval, &
- factor_common,vx,vy,vz,vnspec,is_backward_field)
+ factor_common,vnspec,is_backward_field)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -85,7 +85,7 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: accel_crust_mantle
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
@@ -97,7 +97,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: epsilon_trace_over_3
! [alpha,beta,gamma]val reduced to N_SLS and factor_common to N_SLS*NUM_NODES
- real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
! inner/outer element run flag
@@ -298,7 +298,7 @@
ibool, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
+ one_minus_sum_beta,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc, &
epsilondev_loc, &
@@ -315,7 +315,7 @@
ibool, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
+ one_minus_sum_beta,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc, &
epsilondev_loc, &
@@ -331,7 +331,7 @@
ibool, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
+ one_minus_sum_beta,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc, &
epsilondev_loc, &
@@ -496,7 +496,7 @@
if( ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL ) then
! updates R_memory
call compute_element_att_memory_cm(ispec,R_xx,R_yy,R_xy,R_xz,R_yz, &
- vx,vy,vz,vnspec,factor_common, &
+ vnspec,factor_common, &
alphaval,betaval,gammaval, &
c44store,muvstore, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -36,7 +36,7 @@
epsilondev_xz,epsilondev_yz, &
epsilon_trace_over_3, &
alphaval,betaval,gammaval, &
- factor_common,vx,vy,vz,vnspec)
+ factor_common,vnspec)
use constants_solver
@@ -81,7 +81,7 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: accel_crust_mantle
! variable sized array variables
- integer :: vx,vy,vz,vnspec
+ integer :: vnspec
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
@@ -93,7 +93,7 @@
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: epsilon_trace_over_3
! [alpha,beta,gamma]val reduced to N_SLS and factor_common to N_SLS*NUM_NODES
- real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
! inner/outer element run flag
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -36,7 +36,7 @@
epsilondev_xz,epsilondev_yz, &
epsilon_trace_over_3,&
alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,is_backward_field)
+ vnspec,is_backward_field)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -78,8 +78,8 @@
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
! variable lengths for factor_common and one_minus_sum_beta
- integer vx, vy, vz, vnspec
- real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: factor_common
+ integer :: vnspec
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATT) :: &
@@ -666,7 +666,7 @@
! updates R_memory
call compute_element_att_memory_ic(ispec,R_xx,R_yy,R_xy,R_xz,R_yz, &
- vx,vy,vz,vnspec,factor_common, &
+ vnspec,factor_common, &
alphaval,betaval,gammaval, &
muvstore, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -36,7 +36,7 @@
epsilondev_xz,epsilondev_yz, &
epsilon_trace_over_3,&
alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec)
+ vnspec)
use constants_solver
@@ -79,9 +79,8 @@
! variable lengths for factor_common and one_minus_sum_beta
! variable sized array variables
- integer vx, vy, vz, vnspec
-
- real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+ integer :: vnspec
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATT) :: R_xx,R_yy,R_xy,R_xz,R_yz
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -46,7 +46,7 @@
logical :: phase_is_inner
-!daniel: att - debug
+!daniel debug: att - debug
! integer :: iglob
! logical,parameter :: DEBUG = .false.
! if( DEBUG ) then
@@ -121,9 +121,8 @@
epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
eps_trace_over_3_crust_mantle, &
- 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), .false. )
+ alphaval,betaval,gammaval, &
+ factor_common_crust_mantle,size(factor_common_crust_mantle,5), .false. )
! inner core region
call compute_forces_inner_core_Dev( NSPEC_INNER_CORE_STR_OR_ATT,NGLOB_INNER_CORE, &
@@ -136,9 +135,7 @@
epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
eps_trace_over_3_inner_core,&
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), .false. )
+ factor_common_inner_core,size(factor_common_inner_core,5), .false. )
else
! no Deville optimization
@@ -153,9 +150,8 @@
epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
eps_trace_over_3_crust_mantle, &
- 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) )
+ alphaval,betaval,gammaval, &
+ factor_common_crust_mantle,size(factor_common_crust_mantle,5) )
! inner core region
call compute_forces_inner_core( NSPEC_INNER_CORE_STR_OR_ATT,NGLOB_INNER_CORE, &
NSPEC_INNER_CORE_ATTENUATION, &
@@ -167,9 +163,7 @@
epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
eps_trace_over_3_inner_core,&
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) )
+ factor_common_inner_core,size(factor_common_inner_core,5) )
endif
@@ -189,9 +183,8 @@
b_epsilondev_xy_crust_mantle, &
b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle, &
b_eps_trace_over_3_crust_mantle, &
- 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), .true. )
+ b_alphaval,b_betaval,b_gammaval, &
+ factor_common_crust_mantle,size(factor_common_crust_mantle,5), .true. )
! inner core region
call compute_forces_inner_core_Dev( NSPEC_INNER_CORE_ADJOINT,NGLOB_INNER_CORE_ADJOINT, &
NSPEC_INNER_CORE_STR_AND_ATT, &
@@ -204,9 +197,7 @@
b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core, &
b_eps_trace_over_3_inner_core,&
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), .true. )
+ factor_common_inner_core,size(factor_common_inner_core,5), .true. )
else
! no Deville optimization
@@ -222,9 +213,8 @@
b_epsilondev_xy_crust_mantle, &
b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle, &
b_eps_trace_over_3_crust_mantle, &
- 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) )
+ b_alphaval,b_betaval,b_gammaval, &
+ factor_common_crust_mantle,size(factor_common_crust_mantle,5) )
! inner core region
call compute_forces_inner_core( NSPEC_INNER_CORE_ADJOINT,NGLOB_INNER_CORE_ADJOINT, &
@@ -238,9 +228,7 @@
b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core, &
b_eps_trace_over_3_inner_core,&
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) )
+ factor_common_inner_core,size(factor_common_inner_core,5) )
endif
endif !SIMULATION_TYPE == 3
@@ -526,16 +514,16 @@
! updates (only) acceleration w/ rotation in the crust/mantle region (touches oceans)
if(.NOT. GPU_MODE) then
! on CPU
- call compute_forces_el_update_accel(NGLOB_CRUST_MANTLE,NGLOB_XY_CM,veloc_crust_mantle,accel_crust_mantle, &
- two_omega_earth, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
- NCHUNKS_VAL,ABSORBING_CONDITIONS)
+ call update_accel_elastic(NGLOB_CRUST_MANTLE,NGLOB_XY_CM,veloc_crust_mantle,accel_crust_mantle, &
+ two_omega_earth, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NCHUNKS_VAL,ABSORBING_CONDITIONS)
! adjoint / kernel runs
if (SIMULATION_TYPE == 3) &
- call compute_forces_el_update_accel(NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_XY_CM,b_veloc_crust_mantle,b_accel_crust_mantle, &
- b_two_omega_earth, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
- NCHUNKS_VAL,ABSORBING_CONDITIONS)
+ call update_accel_elastic(NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_XY_CM,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ b_two_omega_earth, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NCHUNKS_VAL,ABSORBING_CONDITIONS)
else
! on GPU
call kernel_3_a_cuda(Mesh_pointer, &
@@ -567,14 +555,14 @@
! (updates velocity)
if(.NOT. GPU_MODE ) then
! on CPU
- call compute_forces_el_update_veloc(NGLOB_CRUST_MANTLE,veloc_crust_mantle,accel_crust_mantle, &
- NGLOB_INNER_CORE,veloc_inner_core,accel_inner_core, &
- deltatover2,two_omega_earth,rmass_inner_core)
+ call update_veloc_elastic(NGLOB_CRUST_MANTLE,veloc_crust_mantle,accel_crust_mantle, &
+ NGLOB_INNER_CORE,veloc_inner_core,accel_inner_core, &
+ deltatover2,two_omega_earth,rmass_inner_core)
! adjoint / kernel runs
if (SIMULATION_TYPE == 3) &
- call compute_forces_el_update_veloc(NGLOB_CRUST_MANTLE_ADJOINT,b_veloc_crust_mantle,b_accel_crust_mantle, &
- NGLOB_INNER_CORE_ADJOINT,b_veloc_inner_core,b_accel_inner_core, &
- b_deltatover2,b_two_omega_earth,rmass_inner_core)
+ call update_veloc_elastic(NGLOB_CRUST_MANTLE_ADJOINT,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ NGLOB_INNER_CORE_ADJOINT,b_veloc_inner_core,b_accel_inner_core, &
+ b_deltatover2,b_two_omega_earth,rmass_inner_core)
else
! on GPU
call kernel_3_b_cuda(Mesh_pointer, &
@@ -582,7 +570,7 @@
endif
-!daniel: att - debug
+!daniel debug: att - debug
! if( DEBUG ) then
! if( SIMULATION_TYPE == 1) then
! if( it > NSTEP - 1000 .and. myrank == 0 ) then
@@ -597,118 +585,3 @@
end subroutine compute_forces_viscoelastic
-
-!=====================================================================
-
- subroutine compute_forces_el_update_accel(NGLOB,NGLOB_XY,veloc_crust_mantle,accel_crust_mantle, &
- two_omega_earth, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
- NCHUNKS_VAL,ABSORBING_CONDITIONS)
-
- use constants_solver,only: CUSTOM_REAL,NDIM
-
- implicit none
-
- integer :: NGLOB,NGLOB_XY,NCHUNKS_VAL
-
- ! velocity & acceleration
- ! crust/mantle region
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc_crust_mantle,accel_crust_mantle
-
- ! mass matrices
- !
- ! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
- ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
- ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
- !
- ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
- ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
- real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmassz_crust_mantle
-
- real(kind=CUSTOM_REAL) :: two_omega_earth
-
- logical :: ABSORBING_CONDITIONS
-
- ! local parameters
- integer :: i
-
- ! updates acceleration w/ rotation in crust/mantle region only
-
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
-
- do i=1,NGLOB
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
- + two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- - 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
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
- + two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
- - two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
- enddo
-
- endif
-
- end subroutine compute_forces_el_update_accel
-
-
-!=====================================================================
-
- subroutine compute_forces_el_update_veloc(NGLOB_CM,veloc_crust_mantle,accel_crust_mantle, &
- NGLOB_IC,veloc_inner_core,accel_inner_core, &
- deltatover2,two_omega_earth,rmass_inner_core)
-
- use constants_solver,only: CUSTOM_REAL,NDIM
-
- implicit none
-
- integer :: NGLOB_CM,NGLOB_IC
-
- ! acceleration & velocity
- ! crust/mantle region
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CM) :: veloc_crust_mantle,accel_crust_mantle
- ! inner core
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_IC) :: veloc_inner_core,accel_inner_core
-
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_IC) :: rmass_inner_core
-
- real(kind=CUSTOM_REAL) :: deltatover2,two_omega_earth
-
- ! local parameters
- integer :: i
-
- ! Newmark time scheme:
- !
- ! note:
- ! - crust/mantle region
- ! needs only velocity corrector terms
- ! (acceleration already updated before)
- ! - inner core region
- ! needs both, acceleration update & velocity corrector terms
-
- ! mantle
- do i=1,NGLOB_CM
- veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
- enddo
- ! inner core
- do i=1,NGLOB_IC
- 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)
-
- veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
- enddo
-
- end subroutine compute_forces_el_update_veloc
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -61,45 +61,87 @@
subroutine compute_kernels_crust_mantle()
use constants_solver
- use specfem_par,only: deltat,GPU_MODE,Mesh_pointer,ANISOTROPIC_KL
+
+ use specfem_par,only: deltat,GPU_MODE,Mesh_pointer,ANISOTROPIC_KL,UNDO_ATTENUATION, &
+ hprime_xx,hprime_xxT,hprime_yy,hprime_zz
+
use specfem_par_crustmantle
+
implicit none
! local parameters
real(kind=CUSTOM_REAL),dimension(21) :: prod
real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: b_epsilondev_loc_matrix
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: b_eps_trace_over_3_loc_matrix
+
integer :: i,j,k,ispec,iglob
if( .not. GPU_MODE ) then
! on CPU
! crust_mantle
do ispec = 1, NSPEC_CRUST_MANTLE
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool_crust_mantle(i,j,k,ispec)
- ! density kernel: see e.g. Tromp et al.(2005), equation (14)
- ! b_displ_crust_mantle is the backward/reconstructed wavefield, that is s(x,t) in eq. (14),
- ! accel_crust_mantle is the adjoint wavefield, that corresponds to s_dagger(x,T-t)
- !
- ! note with respect to eq. (14) the second time derivative is applied to the
- ! adjoint wavefield here rather than the backward/reconstructed wavefield.
- ! this is a valid operation and the resultant kernel identical to the eq. (14).
- !
- ! reason for this is that the adjoint wavefield is in general smoother
- ! since the adjoint sources normally are obtained for filtered traces.
- ! numerically, the time derivative by a finite-difference scheme should
- ! behave better for smoother wavefields, thus containing less numerical artefacts.
- rho_kl_crust_mantle(i,j,k,ispec) = rho_kl_crust_mantle(i,j,k,ispec) &
- + deltat * (accel_crust_mantle(1,iglob) * b_displ_crust_mantle(1,iglob) &
- + accel_crust_mantle(2,iglob) * b_displ_crust_mantle(2,iglob) &
- + accel_crust_mantle(3,iglob) * b_displ_crust_mantle(3,iglob) )
+ ! simulations with UNDO_ATTENUATION save as much memory as possible;
+ ! backward/reconstructed wavefield strain will be re-computed locally here
+ if( UNDO_ATTENUATION ) then
+ if( USE_DEVILLE_PRODUCTS_VAL ) then
+ call compute_element_strain_undo_att_Dev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,&
+ b_displ_crust_mantle,ibool_crust_mantle, &
+ 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,&
+ b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
- ! isotropic kernels
- if (.not. ANISOTROPIC_KL) then
+ else
+ call compute_element_strain_undo_att_noDev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE, &
+ b_displ_crust_mantle, &
+ hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle, &
+ 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, &
+ b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
+ endif
+ else
+ ! backward/reconstructed strain arrays
+ b_eps_trace_over_3_loc_matrix(:,:,:) = b_eps_trace_over_3_crust_mantle(:,:,:,ispec)
+ b_epsilondev_loc_matrix(1,:,:,:) = b_epsilondev_xx_crust_mantle(:,:,:,ispec)
+ b_epsilondev_loc_matrix(2,:,:,:) = b_epsilondev_yy_crust_mantle(:,:,:,ispec)
+ b_epsilondev_loc_matrix(3,:,:,:) = b_epsilondev_xy_crust_mantle(:,:,:,ispec)
+ b_epsilondev_loc_matrix(4,:,:,:) = b_epsilondev_xz_crust_mantle(:,:,:,ispec)
+ b_epsilondev_loc_matrix(5,:,:,:) = b_epsilondev_yz_crust_mantle(:,:,:,ispec)
+ endif
+ ! For anisotropic kernels
+ if (ANISOTROPIC_KL) then
+
+ ! computes fully anisotropic kernel cijkl_kl
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+
+ ! density kernel: see e.g. Tromp et al.(2005), equation (14)
+ ! b_displ_crust_mantle is the backward/reconstructed wavefield, that is s(x,t) in eq. (14),
+ ! accel_crust_mantle is the adjoint wavefield, that corresponds to s_dagger(x,T-t)
+ !
+ ! note with respect to eq. (14) the second time derivative is applied to the
+ ! adjoint wavefield here rather than the backward/reconstructed wavefield.
+ ! this is a valid operation and the resultant kernel identical to the eq. (14).
+ !
+ ! reason for this is that the adjoint wavefield is in general smoother
+ ! since the adjoint sources normally are obtained for filtered traces.
+ ! numerically, the time derivative by a finite-difference scheme should
+ ! behave better for smoother wavefields, thus containing less numerical artefacts.
+ rho_kl_crust_mantle(i,j,k,ispec) = rho_kl_crust_mantle(i,j,k,ispec) &
+ + deltat * (accel_crust_mantle(1,iglob) * b_displ_crust_mantle(1,iglob) &
+ + accel_crust_mantle(2,iglob) * b_displ_crust_mantle(2,iglob) &
+ + accel_crust_mantle(3,iglob) * b_displ_crust_mantle(3,iglob) )
+
+ ! fully anisotropic kernel
! temporary arrays
epsilondev_loc(1) = epsilondev_xx_crust_mantle(i,j,k,ispec)
epsilondev_loc(2) = epsilondev_yy_crust_mantle(i,j,k,ispec)
@@ -107,45 +149,44 @@
epsilondev_loc(4) = epsilondev_xz_crust_mantle(i,j,k,ispec)
epsilondev_loc(5) = epsilondev_yz_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(1) = b_epsilondev_xx_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(2) = b_epsilondev_yy_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(3) = b_epsilondev_xy_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(4) = b_epsilondev_xz_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(5) = b_epsilondev_yz_crust_mantle(i,j,k,ispec)
+ b_epsilondev_loc(:) = b_epsilondev_loc_matrix(:,i,j,k)
- ! kernel for shear modulus, see e.g. Tromp et al. (2005), equation (17)
- ! note: multiplication with 2*mu(x) will be done after the time loop
- beta_kl_crust_mantle(i,j,k,ispec) = beta_kl_crust_mantle(i,j,k,ispec) &
- + deltat * &
- (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
- + (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
- + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
- epsilondev_loc(5)*b_epsilondev_loc(5)) )
+ call compute_strain_product(prod,eps_trace_over_3_crust_mantle(i,j,k,ispec),epsilondev_loc, &
+ b_eps_trace_over_3_loc_matrix(i,j,k),b_epsilondev_loc)
+ cijkl_kl_crust_mantle(:,i,j,k,ispec) = cijkl_kl_crust_mantle(:,i,j,k,ispec) + deltat * prod(:)
- ! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
- ! note: multiplication with kappa(x) will be done after the time loop
- alpha_kl_crust_mantle(i,j,k,ispec) = alpha_kl_crust_mantle(i,j,k,ispec) &
- + deltat * (9 * eps_trace_over_3_crust_mantle(i,j,k,ispec) &
- * b_eps_trace_over_3_crust_mantle(i,j,k,ispec))
-
- endif
-
+ enddo
enddo
enddo
- enddo
- enddo
- ! For anisotropic kernels
- if (ANISOTROPIC_KL) then
+ else
- ! computes fully anisotropic kernel cijkl_kl
- do ispec = 1, NSPEC_CRUST_MANTLE
+ ! isotropic kernels
+
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool_crust_mantle(i,j,k,ispec)
+ ! density kernel: see e.g. Tromp et al.(2005), equation (14)
+ ! b_displ_crust_mantle is the backward/reconstructed wavefield, that is s(x,t) in eq. (14),
+ ! accel_crust_mantle is the adjoint wavefield, that corresponds to s_dagger(x,T-t)
+ !
+ ! note with respect to eq. (14) the second time derivative is applied to the
+ ! adjoint wavefield here rather than the backward/reconstructed wavefield.
+ ! this is a valid operation and the resultant kernel identical to the eq. (14).
+ !
+ ! reason for this is that the adjoint wavefield is in general smoother
+ ! since the adjoint sources normally are obtained for filtered traces.
+ ! numerically, the time derivative by a finite-difference scheme should
+ ! behave better for smoother wavefields, thus containing less numerical artefacts.
+ rho_kl_crust_mantle(i,j,k,ispec) = rho_kl_crust_mantle(i,j,k,ispec) &
+ + deltat * (accel_crust_mantle(1,iglob) * b_displ_crust_mantle(1,iglob) &
+ + accel_crust_mantle(2,iglob) * b_displ_crust_mantle(2,iglob) &
+ + accel_crust_mantle(3,iglob) * b_displ_crust_mantle(3,iglob) )
+
+ ! isotropic kernels
! temporary arrays
epsilondev_loc(1) = epsilondev_xx_crust_mantle(i,j,k,ispec)
epsilondev_loc(2) = epsilondev_yy_crust_mantle(i,j,k,ispec)
@@ -153,24 +194,32 @@
epsilondev_loc(4) = epsilondev_xz_crust_mantle(i,j,k,ispec)
epsilondev_loc(5) = epsilondev_yz_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(1) = b_epsilondev_xx_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(2) = b_epsilondev_yy_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(3) = b_epsilondev_xy_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(4) = b_epsilondev_xz_crust_mantle(i,j,k,ispec)
- b_epsilondev_loc(5) = b_epsilondev_yz_crust_mantle(i,j,k,ispec)
+ b_epsilondev_loc(:) = b_epsilondev_loc_matrix(:,i,j,k)
- call compute_strain_product(prod,eps_trace_over_3_crust_mantle(i,j,k,ispec),epsilondev_loc, &
- b_eps_trace_over_3_crust_mantle(i,j,k,ispec),b_epsilondev_loc)
+ ! kernel for shear modulus, see e.g. Tromp et al. (2005), equation (17)
+ ! note: multiplication with 2*mu(x) will be done after the time loop
+ beta_kl_crust_mantle(i,j,k,ispec) = beta_kl_crust_mantle(i,j,k,ispec) &
+ + deltat * &
+ (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
+ + (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
+ + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
+ epsilondev_loc(5)*b_epsilondev_loc(5)) )
- cijkl_kl_crust_mantle(:,i,j,k,ispec) = cijkl_kl_crust_mantle(:,i,j,k,ispec) + deltat * prod(:)
+ ! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
+ ! note: multiplication with kappa(x) will be done after the time loop
+ alpha_kl_crust_mantle(i,j,k,ispec) = alpha_kl_crust_mantle(i,j,k,ispec) &
+ + deltat * (9 * eps_trace_over_3_crust_mantle(i,j,k,ispec) &
+ * b_eps_trace_over_3_loc_matrix(i,j,k))
+
enddo
enddo
enddo
- enddo
- endif
+ endif ! ANISOTROPIC_KL
+ enddo
+
else
! updates kernel contribution on GPU
@@ -438,8 +487,10 @@
subroutine compute_kernels_inner_core()
use constants_solver
- use specfem_par,only: deltat
- use specfem_par,only: GPU_MODE,Mesh_pointer
+
+ use specfem_par,only: deltat,GPU_MODE,Mesh_pointer,UNDO_ATTENUATION, &
+ hprime_xx,hprime_xxT,hprime_yy,hprime_zz
+
use specfem_par_innercore
implicit none
@@ -448,12 +499,45 @@
real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: b_epsilondev_loc_matrix
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: b_eps_trace_over_3_loc_matrix
+
integer :: i,j,k,ispec,iglob
if( .not. GPU_MODE ) then
! on CPU
! inner_core
do ispec = 1, NSPEC_INNER_CORE
+
+ ! gets element strain
+ if( UNDO_ATTENUATION ) then
+ if( USE_DEVILLE_PRODUCTS_VAL ) then
+ call compute_element_strain_undo_att_Dev(ispec,NGLOB_inner_core,NSPEC_inner_core, &
+ b_displ_inner_core,ibool_inner_core, &
+ 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, &
+ b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
+ else
+ call compute_element_strain_undo_att_noDev(ispec,NGLOB_inner_core,NSPEC_inner_core, &
+ b_displ_inner_core, &
+ hprime_xx,hprime_yy,hprime_zz,ibool_inner_core, &
+ 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, &
+ b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
+ endif
+ else
+ ! backward/reconstructed strain arrays
+ b_eps_trace_over_3_loc_matrix(:,:,:) = b_eps_trace_over_3_inner_core(:,:,:,ispec)
+ b_epsilondev_loc_matrix(1,:,:,:) = b_epsilondev_xx_inner_core(:,:,:,ispec)
+ b_epsilondev_loc_matrix(2,:,:,:) = b_epsilondev_yy_inner_core(:,:,:,ispec)
+ b_epsilondev_loc_matrix(3,:,:,:) = b_epsilondev_xy_inner_core(:,:,:,ispec)
+ b_epsilondev_loc_matrix(4,:,:,:) = b_epsilondev_xz_inner_core(:,:,:,ispec)
+ b_epsilondev_loc_matrix(5,:,:,:) = b_epsilondev_yz_inner_core(:,:,:,ispec)
+ endif
+
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
@@ -470,11 +554,7 @@
epsilondev_loc(4) = epsilondev_xz_inner_core(i,j,k,ispec)
epsilondev_loc(5) = epsilondev_yz_inner_core(i,j,k,ispec)
- b_epsilondev_loc(1) = b_epsilondev_xx_inner_core(i,j,k,ispec)
- b_epsilondev_loc(2) = b_epsilondev_yy_inner_core(i,j,k,ispec)
- b_epsilondev_loc(3) = b_epsilondev_xy_inner_core(i,j,k,ispec)
- b_epsilondev_loc(4) = b_epsilondev_xz_inner_core(i,j,k,ispec)
- b_epsilondev_loc(5) = b_epsilondev_yz_inner_core(i,j,k,ispec)
+ b_epsilondev_loc(:) = b_epsilondev_loc_matrix(:,i,j,k)
beta_kl_inner_core(i,j,k,ispec) = beta_kl_inner_core(i,j,k,ispec) &
+ deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
@@ -485,7 +565,7 @@
alpha_kl_inner_core(i,j,k,ispec) = alpha_kl_inner_core(i,j,k,ispec) &
+ deltat * (9 * eps_trace_over_3_inner_core(i,j,k,ispec) * &
- b_eps_trace_over_3_inner_core(i,j,k,ispec))
+ b_eps_trace_over_3_loc_matrix(i,j,k))
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -208,6 +208,7 @@
NSTEP,it,nit_written)
use constants_solver
+ use specfem_par,only: UNDO_ATTENUATION
implicit none
@@ -270,17 +271,16 @@
double precision, external :: comp_source_time_function
+ ! element strain
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_crust_mantle
+
do irec_local = 1,nrec_local
- ! get global number of that receiver
- irec = number_receiver_global(irec_local)
-
- ! perform the general interpolation using Lagrange polynomials
+ ! initializes
uxd = ZERO
uyd = ZERO
uzd = ZERO
-
-
eps_trace = ZERO
dxx = ZERO
dyy = ZERO
@@ -288,11 +288,38 @@
dxz = ZERO
dyz = ZERO
+ ! gets global number of that receiver
+ irec = number_receiver_global(irec_local)
+
+ ! gets element id
+ ispec = ispec_selected_source(irec)
+
+ ! adjoint strain
+ if( UNDO_ATTENUATION ) then
+ ! recomputes strain
+ call compute_element_strain_undo_att_noDev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE, &
+ displ_crust_mantle, &
+ hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle, &
+ 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, &
+ epsilondev_loc_crust_mantle,eps_trace_over_3_loc_crust_mantle)
+ else
+ ! element adjoint strain
+ eps_trace_over_3_loc_crust_mantle(:,:,:) = eps_trace_over_3_crust_mantle(:,:,:,ispec)
+ epsilondev_loc_crust_mantle(1,:,:,:) = epsilondev_xx_crust_mantle(:,:,:,ispec)
+ epsilondev_loc_crust_mantle(2,:,:,:) = epsilondev_yy_crust_mantle(:,:,:,ispec)
+ epsilondev_loc_crust_mantle(3,:,:,:) = epsilondev_xy_crust_mantle(:,:,:,ispec)
+ epsilondev_loc_crust_mantle(4,:,:,:) = epsilondev_xz_crust_mantle(:,:,:,ispec)
+ epsilondev_loc_crust_mantle(5,:,:,:) = epsilondev_yz_crust_mantle(:,:,:,ispec)
+ endif
+
+ ! perform the general interpolation using Lagrange polynomials
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
- iglob = ibool_crust_mantle(i,j,k,ispec_selected_source(irec))
+ iglob = ibool_crust_mantle(i,j,k,ispec)
hlagrange = hxir_store(irec_local,i)*hetar_store(irec_local,j)*hgammar_store(irec_local,k)
@@ -300,13 +327,14 @@
uyd = uyd + dble(displ_crust_mantle(2,iglob))*hlagrange
uzd = uzd + dble(displ_crust_mantle(3,iglob))*hlagrange
- eps_trace = eps_trace + dble(eps_trace_over_3_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
- dxx = dxx + dble(epsilondev_xx_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
- dyy = dyy + dble(epsilondev_yy_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
- dxy = dxy + dble(epsilondev_xy_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
- dxz = dxz + dble(epsilondev_xz_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
- dyz = dyz + dble(epsilondev_yz_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
+ eps_trace = eps_trace + dble(eps_trace_over_3_loc_crust_mantle(i,j,k))*hlagrange
+ dxx = dxx + dble(epsilondev_loc_crust_mantle(1,i,j,k))*hlagrange
+ dyy = dyy + dble(epsilondev_loc_crust_mantle(2,i,j,k))*hlagrange
+ dxy = dxy + dble(epsilondev_loc_crust_mantle(3,i,j,k))*hlagrange
+ dxz = dxz + dble(epsilondev_loc_crust_mantle(4,i,j,k))*hlagrange
+ dyz = dyz + dble(epsilondev_loc_crust_mantle(5,i,j,k))*hlagrange
+
displ_s(:,i,j,k) = displ_crust_mantle(:,iglob)
enddo
@@ -349,16 +377,14 @@
endif
! frechet derviatives of the source
- ispec = ispec_selected_source(irec)
-
call compute_adj_source_frechet(displ_s,Mxx(irec),Myy(irec),Mzz(irec), &
Mxy(irec),Mxz(irec),Myz(irec),eps_s,eps_m_s,eps_m_l_s, &
hxir_store(irec_local,:),hetar_store(irec_local,:),hgammar_store(irec_local,:), &
hpxir_store(irec_local,:),hpetar_store(irec_local,:),hpgammar_store(irec_local,:), &
hprime_xx,hprime_yy,hprime_zz, &
- xix_crust_mantle(:,:,:,ispec),xiy_crust_mantle(:,:,:,ispec),xiz_crust_mantle(:,:,:,ispec), &
- etax_crust_mantle(:,:,:,ispec),etay_crust_mantle(:,:,:,ispec),etaz_crust_mantle(:,:,:,ispec), &
- gammax_crust_mantle(:,:,:,ispec),gammay_crust_mantle(:,:,:,ispec),gammaz_crust_mantle(:,:,:,ispec))
+ xix_crust_mantle(1,1,1,ispec),xiy_crust_mantle(1,1,1,ispec),xiz_crust_mantle(1,1,1,ispec), &
+ etax_crust_mantle(1,1,1,ispec),etay_crust_mantle(1,1,1,ispec),etaz_crust_mantle(1,1,1,ispec), &
+ gammax_crust_mantle(1,1,1,ispec),gammay_crust_mantle(1,1,1,ispec),gammaz_crust_mantle(1,1,1,ispec))
stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(irec),hdur_gaussian(irec))
stf_deltat = stf * deltat
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -61,7 +61,7 @@
! use the filename to determine the actual contents of the read
if( ADIOS_ENABLED .and. ADIOS_FOR_ARRAYS_SOLVER ) then
call read_attenuation_adios(myrank, prname, &
- factor_common, scale_factor, tau_s, ATT1_VAL, ATT2_VAL, ATT3_VAL, vnspec, T_c_source)
+ factor_common, scale_factor, tau_s, vnspec, T_c_source)
else
open(unit=IIN, file=prname(1:len_trim(prname))//'attenuation.bin', &
status='old',action='read',form='unformatted',iostat=ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -401,6 +401,13 @@
call exit_MPI(myrank, 'error in compiled COMPUTE_AND_STORE_STRAIN parameter, please recompile solver 20')
endif
+ if( FORCE_VECTORIZATION_VAL ) then
+ if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL .and. N_SLS /= 3) &
+ call exit_MPI(myrank, &
+ 'FORCE_VECTORIZATION can only be used with N_SLS == 3 when ATTENUATION .and. .not. PARTIAL_PHYS_DISPERSION_ONLY&
+ & because N_SLS is assumed to be equal to 3 when vectorizing compute_element_iso,tiso,aniso')
+ endif
+
if (SIMULATION_TYPE == 3 .and. (ANISOTROPIC_3D_MANTLE_VAL .or. ANISOTROPIC_INNER_CORE_VAL)) &
call exit_MPI(myrank, 'anisotropic model is not implemented for kernel simulations yet')
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -65,12 +65,6 @@
! get MPI starting time
time_start = wtime()
-#ifdef FORCE_VECTORIZATION
- if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL .and. N_SLS /= 3) &
- stop 'FORCE_VECTORIZATION can only be used with N_SLS == 3 when ATTENUATION .and. .not. PARTIAL_PHYS_DISPERSION_ONLY&
- & because N_SLS is assumed to be equal to 3 when vectorizing compute_element_iso,tiso,aniso'
-#endif
-
! *********************************************************
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
@@ -83,7 +77,7 @@
endif
! update displacement using Newmark time scheme
- call it_update_displacement_scheme()
+ call update_displacement_Newmark()
! acoustic solver for outer core
! (needs to be done first, before elastic one)
@@ -141,128 +135,7 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine it_update_displacement_scheme()
-! explicit Newmark time scheme with acoustic & elastic domains:
-! (see e.g. Hughes, 1987; Chaljub et al., 2003)
-!
-! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
-! chi_dot(t+delta_t) = chi_dot(t) + 1/2 delta_t chi_dot_dot(t) + 1/2 delta_t chi_dot_dot(t+delta_t)
-! chi_dot_dot(t+delta_t) = 1/M_acoustic( -K_acoustic chi(t+delta) + B_acoustic u(t+delta_t) + f(t+delta_t) )
-!
-! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
-! v(t+delta_t) = v(t) + 1/2 delta_t a(t) + 1/2 delta_t a(t+delta_t)
-! a(t+delta_t) = 1/M_elastic ( -K_elastic u(t+delta) + B_elastic chi_dot_dot(t+delta_t) + f( t+delta_t) )
-!
-! where
-! chi, chi_dot, chi_dot_dot are acoustic (fluid) potentials ( dotted with respect to time)
-! u, v, a are displacement,velocity & acceleration
-! M is mass matrix, K stiffness matrix and B boundary term for acoustic/elastic domains
-! f denotes a source term (acoustic/elastic)
-!
-! note that this stage calculates the predictor terms
-!
-! for
-! potential chi_dot(t+delta) requires + 1/2 delta_t chi_dot_dot(t+delta_t)
-! at a later stage (corrector) once where chi_dot_dot(t+delta) is calculated
-! and similar,
-! velocity v(t+delta_t) requires + 1/2 delta_t a(t+delta_t)
-! at a later stage once where a(t+delta) is calculated
-! also:
-! boundary term B_elastic requires chi_dot_dot(t+delta)
-! thus chi_dot_dot has to be updated first before the elastic boundary term is considered
-
- use specfem_par
- use specfem_par_crustmantle
- use specfem_par_innercore
- use specfem_par_outercore
- implicit none
-
- ! local parameters
- integer :: i
-
- ! updates wavefields
- if( .not. GPU_MODE) then
- ! on CPU
-
- ! Newmark time scheme update
- ! mantle
- do i=1,NGLOB_CRUST_MANTLE
- displ_crust_mantle(:,i) = displ_crust_mantle(:,i) &
- + deltat*veloc_crust_mantle(:,i) + deltatsqover2*accel_crust_mantle(:,i)
- veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) &
- + deltatover2*accel_crust_mantle(:,i)
- accel_crust_mantle(:,i) = 0._CUSTOM_REAL
- enddo
- ! outer core
- do i=1,NGLOB_OUTER_CORE
- displ_outer_core(i) = displ_outer_core(i) &
- + deltat*veloc_outer_core(i) + deltatsqover2*accel_outer_core(i)
- veloc_outer_core(i) = veloc_outer_core(i) &
- + deltatover2*accel_outer_core(i)
- accel_outer_core(i) = 0._CUSTOM_REAL
- enddo
- ! inner core
- do i=1,NGLOB_INNER_CORE
- displ_inner_core(:,i) = displ_inner_core(:,i) &
- + deltat*veloc_inner_core(:,i) + deltatsqover2*accel_inner_core(:,i)
- veloc_inner_core(:,i) = veloc_inner_core(:,i) &
- + deltatover2*accel_inner_core(:,i)
- accel_inner_core(:,i) = 0._CUSTOM_REAL
- enddo
-
- ! backward field
- if (SIMULATION_TYPE == 3) then
- ! mantle
- 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)
- b_veloc_crust_mantle(:,i) = b_veloc_crust_mantle(:,i) &
- + b_deltatover2*b_accel_crust_mantle(:,i)
- b_accel_crust_mantle(:,i) = 0._CUSTOM_REAL
- enddo
- ! outer core
- do i=1,NGLOB_OUTER_CORE
- b_displ_outer_core(i) = b_displ_outer_core(i) &
- + b_deltat*b_veloc_outer_core(i) + b_deltatsqover2*b_accel_outer_core(i)
- b_veloc_outer_core(i) = b_veloc_outer_core(i) &
- + b_deltatover2*b_accel_outer_core(i)
- b_accel_outer_core(i) = 0._CUSTOM_REAL
- enddo
- ! inner core
- 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)
- b_veloc_inner_core(:,i) = b_veloc_inner_core(:,i) &
- + b_deltatover2*b_accel_inner_core(:,i)
- b_accel_inner_core(:,i) = 0._CUSTOM_REAL
- enddo
- endif ! SIMULATION_TYPE == 3
- else
- ! on GPU
- ! Includes SIM_TYPE 1 & 3
-
- ! outer core region
- call it_update_displacement_oc_cuda(Mesh_pointer, &
- deltat, deltatsqover2, deltatover2, &
- b_deltat, b_deltatsqover2, b_deltatover2)
- ! inner core region
- call it_update_displacement_ic_cuda(Mesh_pointer, &
- deltat, deltatsqover2, deltatover2, &
- b_deltat, b_deltatsqover2, b_deltatover2)
-
- ! crust/mantle region
- call it_update_displacement_cm_cuda(Mesh_pointer, &
- deltat, deltatsqover2, deltatover2, &
- b_deltat, b_deltatsqover2, b_deltatover2)
- endif
-
- end subroutine it_update_displacement_scheme
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
subroutine it_print_elapsed_time()
use specfem_par
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -557,6 +557,12 @@
call gather_all_i(ispec_selected_source_subset,NSOURCES_SUBSET_current_size, &
ispec_selected_source_all,NSOURCES_SUBSET_current_size,NPROCTOT_VAL)
+ ! checks that the gather operation went well
+ if(myrank == 0) then
+ if(minval(ispec_selected_source_all(:,:)) <= 0) &
+ call exit_MPI(myrank,'gather operation failed for source')
+ endif
+
call gather_all_dp(xi_source_subset,NSOURCES_SUBSET_current_size, &
xi_source_all,NSOURCES_SUBSET_current_size,NPROCTOT_VAL)
call gather_all_dp(eta_source_subset,NSOURCES_SUBSET_current_size, &
@@ -575,10 +581,6 @@
! this is executed by main process only
if(myrank == 0) then
- ! check that the gather operation went well
- if(minval(ispec_selected_source_all) <= 0) &
- call exit_MPI(myrank,'gather operation failed for source')
-
! loop on all the sources within subsets
do isource_in_this_subset = 1,NSOURCES_SUBSET_current_size
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -183,8 +183,10 @@
write(IMAIN,*)
if(ATTENUATION_VAL) then
write(IMAIN,*) 'incorporating attenuation using ',N_SLS,' standard linear solids'
- if(ATTENUATION_3D_VAL) write(IMAIN,*) 'using 3D attenuation model'
- if(PARTIAL_PHYS_DISPERSION_ONLY ) write(IMAIN,*) 'mimicking effects on velocity only'
+ if(ATTENUATION_3D_VAL) &
+ write(IMAIN,*) 'using 3D attenuation model'
+ if(PARTIAL_PHYS_DISPERSION_ONLY_VAL ) &
+ write(IMAIN,*) 'mimicking effects on velocity only'
else
write(IMAIN,*) 'no attenuation'
endif
@@ -564,7 +566,7 @@
call flush_IMAIN()
endif
-
+ ! checks
! the following has to be true for the the array dimensions of eps to match with those of xstore etc..
! note that epsilondev and eps_trace_over_3 don't have the same dimensions.. could cause trouble
if (NSPEC_CRUST_MANTLE_STR_OR_ATT /= NSPEC_CRUST_MANTLE) &
@@ -572,6 +574,7 @@
if (NSPEC_CRUST_MANTLE_STRAIN_ONLY /= NSPEC_CRUST_MANTLE) &
stop 'NSPEC_CRUST_MANTLE_STRAIN_ONLY /= NSPEC_CRUST_MANTLE'
+ ! counts total number of points for movie file output
call movie_volume_count_points()
allocate(nu_3dmovie(3,3,npoints_3dmovie),stat=ier)
@@ -807,7 +810,8 @@
use specfem_par
use specfem_par_crustmantle
use specfem_par_innercore
- use specfem_par_movie
+ use specfem_par_movie,only: muvstore_crust_mantle_3dmovie
+
implicit none
! local parameters
@@ -903,8 +907,8 @@
+ scale_factor_minus_one * mul
else
if(MOVIE_VOLUME .and. SIMULATION_TYPE==3) then
- ! store the original value of \mu to comput \mu*\eps
- muvstore_crust_mantle_3dmovie(i,j,k,ispec)=muvstore_crust_mantle(i,j,k,ispec)
+ ! store the original value of \mu to compute \mu*\eps
+ muvstore_crust_mantle_3dmovie(i,j,k,ispec) = muvstore_crust_mantle(i,j,k,ispec)
endif
muvstore_crust_mantle(i,j,k,ispec) = muvstore_crust_mantle(i,j,k,ispec) * scale_factor
@@ -1085,39 +1089,75 @@
endif
! initialize to be on the save side for adjoint runs SIMULATION_TYPE==2
+ ! crust/mantle
eps_trace_over_3_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
-
epsilondev_xx_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yy_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xy_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xz_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yz_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
+ ! backward/reconstructed strain fields
+ if( .not. UNDO_ATTENUATION ) then
+ allocate(b_epsilondev_xx_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ b_epsilondev_yy_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ b_epsilondev_xy_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ b_epsilondev_xz_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ b_epsilondev_yz_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ b_eps_trace_over_3_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_epsilondev*** arrays for crust/mantle')
+
+ allocate(b_epsilondev_xx_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT), &
+ b_epsilondev_yy_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT), &
+ b_epsilondev_xy_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT), &
+ b_epsilondev_xz_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT), &
+ b_epsilondev_yz_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT), &
+ b_eps_trace_over_3_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_epsilondev*** arrays for inner core')
+ else
+ ! dummy arrays
+ allocate(b_epsilondev_xx_crust_mantle(1,1,1,1), &
+ b_epsilondev_yy_crust_mantle(1,1,1,1), &
+ b_epsilondev_xy_crust_mantle(1,1,1,1), &
+ b_epsilondev_xz_crust_mantle(1,1,1,1), &
+ b_epsilondev_yz_crust_mantle(1,1,1,1), &
+ b_eps_trace_over_3_crust_mantle(1,1,1,1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_epsilondev*** arrays for crust/mantle')
+
+ allocate(b_epsilondev_xx_inner_core(1,1,1,1), &
+ b_epsilondev_yy_inner_core(1,1,1,1), &
+ b_epsilondev_xy_inner_core(1,1,1,1), &
+ b_epsilondev_xz_inner_core(1,1,1,1), &
+ b_epsilondev_yz_inner_core(1,1,1,1), &
+ b_eps_trace_over_3_inner_core(1,1,1,1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_epsilondev*** arrays for inner core')
+ endif
+
+ ! inner core
eps_trace_over_3_inner_core(:,:,:,:) = 0._CUSTOM_REAL
-
epsilondev_xx_inner_core(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yy_inner_core(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xy_inner_core(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xz_inner_core(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yz_inner_core(:,:,:,:) = 0._CUSTOM_REAL
+
+
if(FIX_UNDERFLOW_PROBLEM) then
+ ! crust/mantle
eps_trace_over_3_crust_mantle(:,:,:,:) = VERYSMALLVAL
-
epsilondev_xx_crust_mantle(:,:,:,:) = VERYSMALLVAL
epsilondev_yy_crust_mantle(:,:,:,:) = VERYSMALLVAL
epsilondev_xy_crust_mantle(:,:,:,:) = VERYSMALLVAL
epsilondev_xz_crust_mantle(:,:,:,:) = VERYSMALLVAL
epsilondev_yz_crust_mantle(:,:,:,:) = VERYSMALLVAL
-
+ ! inner core
eps_trace_over_3_inner_core(:,:,:,:) = VERYSMALLVAL
-
epsilondev_xx_inner_core(:,:,:,:) = VERYSMALLVAL
epsilondev_yy_inner_core(:,:,:,:) = VERYSMALLVAL
epsilondev_xy_inner_core(:,:,:,:) = VERYSMALLVAL
epsilondev_xz_inner_core(:,:,:,:) = VERYSMALLVAL
epsilondev_yz_inner_core(:,:,:,:) = VERYSMALLVAL
-
endif
if (COMPUTE_AND_STORE_STRAIN) then
@@ -1556,7 +1596,7 @@
OCEANS_VAL, &
GRAVITY_VAL, &
ROTATION_VAL, &
- ATTENUATION_VAL, &
+ ATTENUATION_VAL,UNDO_ATTENUATION, &
PARTIAL_PHYS_DISPERSION_ONLY,USE_3D_ATTENUATION_ARRAYS, &
COMPUTE_AND_STORE_STRAIN, &
ANISOTROPIC_3D_MANTLE_VAL,ANISOTROPIC_INNER_CORE_VAL, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_attenuation_adios.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_attenuation_adios.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_attenuation_adios.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -30,10 +30,9 @@
!> \brief Read adios attenuation arrays created by the mesher
! (regX_attenuation.bp)
subroutine read_attenuation_adios(myrank, prname, &
- factor_common, scale_factor, tau_s, vx, vy, vz, vnspec, T_c_source)
+ factor_common, scale_factor, tau_s, vnspec, T_c_source)
- use constants
-
+ use constants_solver
use mpi
use adios_read_mod
use specfem_par,only: ATTENUATION_VAL
@@ -42,9 +41,9 @@
integer :: myrank
- integer :: vx,vy,vz,vnspec
- real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: scale_factor
- real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+ integer :: vnspec
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: scale_factor
+ real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: factor_common
double precision, dimension(N_SLS) :: tau_s
character(len=150) :: prname
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -37,85 +37,79 @@
implicit none
! local parameters
+ integer :: ier
character(len=150) outputname
+ ! checks
+ if(NUMBER_OF_RUNS < 1 .or. NUMBER_OF_RUNS > NSTEP) &
+ stop 'number of restart runs can not be less than 1 or greater than NSTEP'
+ if(NUMBER_OF_THIS_RUN < 1 .or. NUMBER_OF_THIS_RUN > NUMBER_OF_RUNS) &
+ stop 'incorrect run number'
+ if (SIMULATION_TYPE /= 1 .and. NUMBER_OF_RUNS /= 1) &
+ stop 'Only 1 run for SIMULATION_TYPE = 2/3'
+
! define correct time steps if restart files
- if(NUMBER_OF_RUNS < 1 .or. NUMBER_OF_RUNS > 3) stop 'number of restart runs can be 1, 2 or 3'
- if(NUMBER_OF_THIS_RUN < 1 .or. NUMBER_OF_THIS_RUN > NUMBER_OF_RUNS) stop 'incorrect run number'
- if (SIMULATION_TYPE /= 1 .and. NUMBER_OF_RUNS /= 1) stop 'Only 1 run for SIMULATION_TYPE = 2/3'
-
! set start/end steps for time iteration loop
- select case( NUMBER_OF_RUNS )
- case( 3 )
- if(NUMBER_OF_THIS_RUN == 1) then
- it_begin = 1
- it_end = NSTEP/3
- else if(NUMBER_OF_THIS_RUN == 2) then
- it_begin = NSTEP/3 + 1
- it_end = 2*(NSTEP/3)
- else
- it_begin = 2*(NSTEP/3) + 1
- it_end = NSTEP
- endif
- case( 2 )
- if(NUMBER_OF_THIS_RUN == 1) then
- it_begin = 1
- it_end = NSTEP/2
- else
- it_begin = NSTEP/2 + 1
- it_end = NSTEP
- endif
- case default
- it_begin = 1
+ it_begin = (NUMBER_OF_THIS_RUN - 1) * (NSTEP / NUMBER_OF_RUNS) + 1
+ if (NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS) then
+ it_end = NUMBER_OF_THIS_RUN * (NSTEP / NUMBER_OF_RUNS)
+ else
+ ! Last run may be a bit larger
it_end = NSTEP
- end select
+ endif
+ ! checks if anything to do
+ if( UNDO_ATTENUATION ) return
+
! read files back from local disk or MT tape system if restart file
if(NUMBER_OF_THIS_RUN > 1) then
if( ADIOS_ENABLED .and. ADIOS_FOR_FORWARD_ARRAYS ) then
call read_intermediate_forward_arrays_adios()
else
write(outputname,"('dump_all_arrays',i6.6)") myrank
- open(unit=55,file=trim(LOCAL_TMP_PATH)//'/'//outputname,status='old',action='read',form='unformatted')
+ open(unit=IIN,file=trim(LOCAL_TMP_PATH)//'/'//outputname,status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file dump_all_arrays*** for reading')
- read(55) displ_crust_mantle
- read(55) veloc_crust_mantle
- read(55) accel_crust_mantle
- read(55) displ_inner_core
- read(55) veloc_inner_core
- read(55) accel_inner_core
- read(55) displ_outer_core
- read(55) veloc_outer_core
- read(55) accel_outer_core
+ read(IIN) displ_crust_mantle
+ read(IIN) veloc_crust_mantle
+ read(IIN) accel_crust_mantle
- read(55) epsilondev_xx_crust_mantle
- read(55) epsilondev_yy_crust_mantle
- read(55) epsilondev_xy_crust_mantle
- read(55) epsilondev_xz_crust_mantle
- read(55) epsilondev_yz_crust_mantle
+ read(IIN) displ_inner_core
+ read(IIN) veloc_inner_core
+ read(IIN) accel_inner_core
- read(55) epsilondev_xx_inner_core
- read(55) epsilondev_yy_inner_core
- read(55) epsilondev_xy_inner_core
- read(55) epsilondev_xz_inner_core
- read(55) epsilondev_yz_inner_core
+ read(IIN) displ_outer_core
+ read(IIN) veloc_outer_core
+ read(IIN) accel_outer_core
- read(55) A_array_rotation
- read(55) B_array_rotation
+ read(IIN) epsilondev_xx_crust_mantle
+ read(IIN) epsilondev_yy_crust_mantle
+ read(IIN) epsilondev_xy_crust_mantle
+ read(IIN) epsilondev_xz_crust_mantle
+ read(IIN) epsilondev_yz_crust_mantle
- read(55) R_xx_crust_mantle
- read(55) R_yy_crust_mantle
- read(55) R_xy_crust_mantle
- read(55) R_xz_crust_mantle
- read(55) R_yz_crust_mantle
+ read(IIN) epsilondev_xx_inner_core
+ read(IIN) epsilondev_yy_inner_core
+ read(IIN) epsilondev_xy_inner_core
+ read(IIN) epsilondev_xz_inner_core
+ read(IIN) epsilondev_yz_inner_core
- read(55) R_xx_inner_core
- read(55) R_yy_inner_core
- read(55) R_xy_inner_core
- read(55) R_xz_inner_core
- read(55) R_yz_inner_core
+ read(IIN) A_array_rotation
+ read(IIN) B_array_rotation
- close(55)
+ read(IIN) R_xx_crust_mantle
+ read(IIN) R_yy_crust_mantle
+ read(IIN) R_xy_crust_mantle
+ read(IIN) R_xz_crust_mantle
+ read(IIN) R_yz_crust_mantle
+
+ read(IIN) R_xx_inner_core
+ read(IIN) R_yy_inner_core
+ read(IIN) R_xy_inner_core
+ read(IIN) R_xz_inner_core
+ read(IIN) R_yz_inner_core
+
+ close(IIN)
endif
endif
@@ -191,7 +185,7 @@
call read_forward_arrays_adios()
else
write(outputname,'(a,i6.6,a)') 'proc',myrank,'_save_forward_arrays.bin'
- open(unit=55,file=trim(LOCAL_TMP_PATH)//'/'//outputname, &
+ open(unit=IIN,file=trim(LOCAL_TMP_PATH)//'/'//outputname, &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error: opening proc_****_save_forward_arrays.bin'
@@ -199,49 +193,51 @@
call exit_mpi(myrank,'error open file save_forward_arrays.bin')
endif
- read(55) b_displ_crust_mantle
- read(55) b_veloc_crust_mantle
- read(55) b_accel_crust_mantle
+ read(IIN) b_displ_crust_mantle
+ read(IIN) b_veloc_crust_mantle
+ read(IIN) b_accel_crust_mantle
- read(55) b_displ_inner_core
- read(55) b_veloc_inner_core
- read(55) b_accel_inner_core
+ read(IIN) b_displ_inner_core
+ read(IIN) b_veloc_inner_core
+ read(IIN) b_accel_inner_core
- read(55) b_displ_outer_core
- read(55) b_veloc_outer_core
- read(55) b_accel_outer_core
+ read(IIN) b_displ_outer_core
+ read(IIN) b_veloc_outer_core
+ read(IIN) b_accel_outer_core
- read(55) b_epsilondev_xx_crust_mantle
- read(55) b_epsilondev_yy_crust_mantle
- read(55) b_epsilondev_xy_crust_mantle
- read(55) b_epsilondev_xz_crust_mantle
- read(55) b_epsilondev_yz_crust_mantle
+ if( .not. UNDO_ATTENUATION ) then
+ read(IIN) b_epsilondev_xx_crust_mantle
+ read(IIN) b_epsilondev_yy_crust_mantle
+ read(IIN) b_epsilondev_xy_crust_mantle
+ read(IIN) b_epsilondev_xz_crust_mantle
+ read(IIN) b_epsilondev_yz_crust_mantle
- read(55) b_epsilondev_xx_inner_core
- read(55) b_epsilondev_yy_inner_core
- read(55) b_epsilondev_xy_inner_core
- read(55) b_epsilondev_xz_inner_core
- read(55) b_epsilondev_yz_inner_core
+ read(IIN) b_epsilondev_xx_inner_core
+ read(IIN) b_epsilondev_yy_inner_core
+ read(IIN) b_epsilondev_xy_inner_core
+ read(IIN) b_epsilondev_xz_inner_core
+ read(IIN) b_epsilondev_yz_inner_core
+ endif
if (ROTATION_VAL) then
- read(55) b_A_array_rotation
- read(55) b_B_array_rotation
+ read(IIN) b_A_array_rotation
+ read(IIN) b_B_array_rotation
endif
if (ATTENUATION_VAL) then
- read(55) b_R_xx_crust_mantle
- read(55) b_R_yy_crust_mantle
- read(55) b_R_xy_crust_mantle
- read(55) b_R_xz_crust_mantle
- read(55) b_R_yz_crust_mantle
+ read(IIN) b_R_xx_crust_mantle
+ read(IIN) b_R_yy_crust_mantle
+ read(IIN) b_R_xy_crust_mantle
+ read(IIN) b_R_xz_crust_mantle
+ read(IIN) b_R_yz_crust_mantle
- read(55) b_R_xx_inner_core
- read(55) b_R_yy_inner_core
- read(55) b_R_xy_inner_core
- read(55) b_R_xz_inner_core
- read(55) b_R_yz_inner_core
+ read(IIN) b_R_xx_inner_core
+ read(IIN) b_R_yy_inner_core
+ read(IIN) b_R_xy_inner_core
+ read(IIN) b_R_xz_inner_core
+ read(IIN) b_R_yz_inner_core
endif
- close(55)
+ close(IIN)
endif ! ADIOS_FOR_FORWARD_ARRAYS
! transfers fields onto GPU
@@ -258,15 +254,18 @@
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
Mesh_pointer)
- call transfer_b_strain_cm_to_device(Mesh_pointer, &
+ if( .not. UNDO_ATTENUATION ) then
+ call transfer_b_strain_cm_to_device(Mesh_pointer, &
b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle, &
b_epsilondev_xy_crust_mantle,b_epsilondev_xz_crust_mantle, &
b_epsilondev_yz_crust_mantle)
- call transfer_b_strain_ic_to_device(Mesh_pointer, &
+ call transfer_b_strain_ic_to_device(Mesh_pointer, &
b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core, &
b_epsilondev_xy_inner_core,b_epsilondev_xz_inner_core, &
b_epsilondev_yz_inner_core)
+ endif
+
if (ROTATION_VAL) then
call transfer_b_rotation_to_device(Mesh_pointer,b_A_array_rotation,b_B_array_rotation)
endif
@@ -284,3 +283,64 @@
endif
end subroutine read_forward_arrays
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine read_forward_arrays_undoatt(iteration_on_subset)
+
+! reads in saved wavefields
+
+ use specfem_par
+ use specfem_par_crustmantle
+ use specfem_par_innercore
+ use specfem_par_outercore
+
+ implicit none
+
+ integer :: iteration_on_subset
+
+ ! local parameters
+ integer :: ier
+ character(len=150) :: outputname
+
+ write(outputname,'(a,i6.6,a,i6.6,a)') 'proc',myrank,'_save_frame_at',iteration_on_subset,'.bin'
+ open(unit=IIN,file=trim(LOCAL_PATH)//'/'//outputname,status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file proc***_save_frame_at** for reading')
+
+ read(IIN) b_displ_crust_mantle
+ read(IIN) b_veloc_crust_mantle
+ read(IIN) b_accel_crust_mantle
+
+ read(IIN) b_displ_inner_core
+ read(IIN) b_veloc_inner_core
+ read(IIN) b_accel_inner_core
+
+ read(IIN) b_displ_outer_core
+ read(IIN) b_veloc_outer_core
+ read(IIN) b_accel_outer_core
+
+ if (ROTATION_VAL) then
+ read(IIN) b_A_array_rotation
+ read(IIN) b_B_array_rotation
+ endif
+
+ if (ATTENUATION_VAL) then
+ read(IIN) b_R_xx_crust_mantle
+ read(IIN) b_R_yy_crust_mantle
+ read(IIN) b_R_xy_crust_mantle
+ read(IIN) b_R_xz_crust_mantle
+ read(IIN) b_R_yz_crust_mantle
+
+ read(IIN) b_R_xx_inner_core
+ read(IIN) b_R_yy_inner_core
+ read(IIN) b_R_xy_inner_core
+ read(IIN) b_R_xz_inner_core
+ read(IIN) b_R_yz_inner_core
+ endif
+
+ close(IIN)
+
+ end subroutine read_forward_arrays_undoatt
+
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-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk 2013-08-29 07:47:36 UTC (rev 22735)
@@ -59,6 +59,7 @@
$O/compute_boundary_kernel.solverstatic.o \
$O/compute_coupling.solverstatic.o \
$O/compute_element.solverstatic.o \
+ $O/compute_element_strain.solverstatic.o \
$O/compute_forces_acoustic_calling_routine.solverstatic.o \
$O/compute_forces_viscoelastic_calling_routine.solverstatic.o \
$O/compute_forces_crust_mantle_noDev.solverstatic.o \
@@ -90,6 +91,7 @@
$O/setup_GLL_points.solverstatic.o \
$O/setup_sources_receivers.solverstatic.o \
$O/specfem3D.solverstatic.o \
+ $O/update_displacement_Newmark.solverstatic.o \
$O/write_movie_output.solverstatic.o \
$O/write_movie_volume.solverstatic.o \
$O/write_movie_surface.solverstatic.o \
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -34,7 +34,7 @@
!! Moreover, one can choose to save these arrays with the help of the ADIOS
!! library. Set ADIOS_FOR_FORWARD_ARRAYS to true in the Par_file for this
!! feature.
-subroutine save_forward_arrays()
+ subroutine save_forward_arrays()
use specfem_par
use specfem_par_crustmantle
@@ -44,55 +44,62 @@
implicit none
! local parameters
+ integer :: ier
character(len=150) outputname
+ ! checks if anything to do
+ if( UNDO_ATTENUATION ) return
+
! save files to local disk or tape system if restart file
if(NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS) then
if( ADIOS_ENABLED .and. ADIOS_FOR_FORWARD_ARRAYS ) then
call save_intermediate_forward_arrays_adios()
else
write(outputname,"('dump_all_arrays',i6.6)") myrank
- open(unit=55,file=trim(LOCAL_TMP_PATH)//'/'//outputname, &
- status='unknown',form='unformatted',action='write')
+ open(unit=IOUT,file=trim(LOCAL_TMP_PATH)//'/'//outputname, &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file dump_all_arrays*** for writing')
- write(55) displ_crust_mantle
- write(55) veloc_crust_mantle
- write(55) accel_crust_mantle
- write(55) displ_inner_core
- write(55) veloc_inner_core
- write(55) accel_inner_core
- write(55) displ_outer_core
- write(55) veloc_outer_core
- write(55) accel_outer_core
+ write(IOUT) displ_crust_mantle
+ write(IOUT) veloc_crust_mantle
+ write(IOUT) accel_crust_mantle
- write(55) epsilondev_xx_crust_mantle
- write(55) epsilondev_yy_crust_mantle
- write(55) epsilondev_xy_crust_mantle
- write(55) epsilondev_xz_crust_mantle
- write(55) epsilondev_yz_crust_mantle
+ write(IOUT) displ_inner_core
+ write(IOUT) veloc_inner_core
+ write(IOUT) accel_inner_core
- write(55) epsilondev_xx_inner_core
- write(55) epsilondev_yy_inner_core
- write(55) epsilondev_xy_inner_core
- write(55) epsilondev_xz_inner_core
- write(55) epsilondev_yz_inner_core
+ write(IOUT) displ_outer_core
+ write(IOUT) veloc_outer_core
+ write(IOUT) accel_outer_core
- write(55) A_array_rotation
- write(55) B_array_rotation
+ write(IOUT) epsilondev_xx_crust_mantle
+ write(IOUT) epsilondev_yy_crust_mantle
+ write(IOUT) epsilondev_xy_crust_mantle
+ write(IOUT) epsilondev_xz_crust_mantle
+ write(IOUT) epsilondev_yz_crust_mantle
- write(55) R_xx_crust_mantle
- write(55) R_yy_crust_mantle
- write(55) R_xy_crust_mantle
- write(55) R_xz_crust_mantle
- write(55) R_yz_crust_mantle
+ write(IOUT) epsilondev_xx_inner_core
+ write(IOUT) epsilondev_yy_inner_core
+ write(IOUT) epsilondev_xy_inner_core
+ write(IOUT) epsilondev_xz_inner_core
+ write(IOUT) epsilondev_yz_inner_core
- write(55) R_xx_inner_core
- write(55) R_yy_inner_core
- write(55) R_xy_inner_core
- write(55) R_xz_inner_core
- write(55) R_yz_inner_core
+ write(IOUT) A_array_rotation
+ write(IOUT) B_array_rotation
- close(55)
+ write(IOUT) R_xx_crust_mantle
+ write(IOUT) R_yy_crust_mantle
+ write(IOUT) R_xy_crust_mantle
+ write(IOUT) R_xz_crust_mantle
+ write(IOUT) R_yz_crust_mantle
+
+ write(IOUT) R_xx_inner_core
+ write(IOUT) R_yy_inner_core
+ write(IOUT) R_xy_inner_core
+ write(IOUT) R_xz_inner_core
+ write(IOUT) R_yz_inner_core
+
+ close(IOUT)
endif
endif
@@ -102,54 +109,124 @@
call save_forward_arrays_adios()
else
write(outputname,'(a,i6.6,a)') 'proc',myrank,'_save_forward_arrays.bin'
- open(unit=55,file=trim(LOCAL_TMP_PATH)//'/'//outputname,status='unknown', &
- form='unformatted',action='write')
+ open(unit=IOUT,file=trim(LOCAL_TMP_PATH)//'/'//outputname,status='unknown', &
+ form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file proc***_save_forward_arrays** for writing')
- write(55) displ_crust_mantle
- write(55) veloc_crust_mantle
- write(55) accel_crust_mantle
+ write(IOUT) displ_crust_mantle
+ write(IOUT) veloc_crust_mantle
+ write(IOUT) accel_crust_mantle
- write(55) displ_inner_core
- write(55) veloc_inner_core
- write(55) accel_inner_core
+ write(IOUT) displ_inner_core
+ write(IOUT) veloc_inner_core
+ write(IOUT) accel_inner_core
- write(55) displ_outer_core
- write(55) veloc_outer_core
- write(55) accel_outer_core
+ write(IOUT) displ_outer_core
+ write(IOUT) veloc_outer_core
+ write(IOUT) accel_outer_core
- write(55) epsilondev_xx_crust_mantle
- write(55) epsilondev_yy_crust_mantle
- write(55) epsilondev_xy_crust_mantle
- write(55) epsilondev_xz_crust_mantle
- write(55) epsilondev_yz_crust_mantle
+ if( .not. UNDO_ATTENUATION ) then
+ write(IOUT) epsilondev_xx_crust_mantle
+ write(IOUT) epsilondev_yy_crust_mantle
+ write(IOUT) epsilondev_xy_crust_mantle
+ write(IOUT) epsilondev_xz_crust_mantle
+ write(IOUT) epsilondev_yz_crust_mantle
- write(55) epsilondev_xx_inner_core
- write(55) epsilondev_yy_inner_core
- write(55) epsilondev_xy_inner_core
- write(55) epsilondev_xz_inner_core
- write(55) epsilondev_yz_inner_core
+ write(IOUT) epsilondev_xx_inner_core
+ write(IOUT) epsilondev_yy_inner_core
+ write(IOUT) epsilondev_xy_inner_core
+ write(IOUT) epsilondev_xz_inner_core
+ write(IOUT) epsilondev_yz_inner_core
+ endif
if (ROTATION_VAL) then
- write(55) A_array_rotation
- write(55) B_array_rotation
+ write(IOUT) A_array_rotation
+ write(IOUT) B_array_rotation
endif
if (ATTENUATION_VAL) then
- write(55) R_xx_crust_mantle
- write(55) R_yy_crust_mantle
- write(55) R_xy_crust_mantle
- write(55) R_xz_crust_mantle
- write(55) R_yz_crust_mantle
+ write(IOUT) R_xx_crust_mantle
+ write(IOUT) R_yy_crust_mantle
+ write(IOUT) R_xy_crust_mantle
+ write(IOUT) R_xz_crust_mantle
+ write(IOUT) R_yz_crust_mantle
- write(55) R_xx_inner_core
- write(55) R_yy_inner_core
- write(55) R_xy_inner_core
- write(55) R_xz_inner_core
- write(55) R_yz_inner_core
+ write(IOUT) R_xx_inner_core
+ write(IOUT) R_yy_inner_core
+ write(IOUT) R_xy_inner_core
+ write(IOUT) R_xz_inner_core
+ write(IOUT) R_yz_inner_core
endif
- close(55)
+ close(IOUT)
endif
endif
-end subroutine save_forward_arrays
+ end subroutine save_forward_arrays
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine save_forward_arrays_undoatt(iteration_on_subset)
+
+ use specfem_par
+ use specfem_par_crustmantle
+ use specfem_par_innercore
+ use specfem_par_outercore
+
+ implicit none
+
+ integer :: iteration_on_subset
+
+ ! local parameters
+ integer :: ier
+ character(len=150) :: outputname
+
+ ! save files to local disk or tape system if restart file
+ if(NUMBER_OF_RUNS > 1) stop 'NUMBER_OF_RUNS > 1 is not support for undoing attenuation'
+
+ ! save last frame of the forward simulation
+ if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
+ write(outputname,'(a,i6.6,a,i6.6,a)') 'proc',myrank,'_save_frame_at',iteration_on_subset,'.bin'
+ open(unit=IOUT,file=trim(LOCAL_PATH)//'/'//outputname,status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file proc***_save_frame_at** for writing')
+
+ write(IOUT) displ_crust_mantle
+ write(IOUT) veloc_crust_mantle
+ write(IOUT) accel_crust_mantle
+
+ write(IOUT) displ_inner_core
+ write(IOUT) veloc_inner_core
+ write(IOUT) accel_inner_core
+
+ write(IOUT) displ_outer_core
+ write(IOUT) veloc_outer_core
+ write(IOUT) accel_outer_core
+
+ if (ROTATION_VAL) then
+ write(IOUT) A_array_rotation
+ write(IOUT) B_array_rotation
+ endif
+
+ if (ATTENUATION_VAL) then
+ write(IOUT) R_xx_crust_mantle
+ write(IOUT) R_yy_crust_mantle
+ write(IOUT) R_xy_crust_mantle
+ write(IOUT) R_xz_crust_mantle
+ write(IOUT) R_yz_crust_mantle
+
+ write(IOUT) R_xx_inner_core
+ write(IOUT) R_yy_inner_core
+ write(IOUT) R_xy_inner_core
+ write(IOUT) R_xz_inner_core
+ write(IOUT) R_yz_inner_core
+ endif
+
+ close(IOUT)
+
+ endif
+
+ end subroutine save_forward_arrays_undoatt
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-08-28 17:12:17 UTC (rev 22734)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -369,13 +369,21 @@
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_AND_ATT) :: &
b_R_xx_crust_mantle,b_R_yy_crust_mantle,b_R_xy_crust_mantle,b_R_xz_crust_mantle,b_R_yz_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle,b_epsilondev_xy_crust_mantle, &
b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
b_eps_trace_over_3_crust_mantle
+! daniel debug: static
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+! b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle,b_epsilondev_xy_crust_mantle, &
+! b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle
+
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+! b_eps_trace_over_3_crust_mantle
+
! for crust/oceans coupling
integer, dimension(NSPEC2DMAX_XMIN_XMAX_CM) :: ibelm_xmin_crust_mantle,ibelm_xmax_crust_mantle
integer, dimension(NSPEC2DMAX_YMIN_YMAX_CM) :: ibelm_ymin_crust_mantle,ibelm_ymax_crust_mantle
@@ -526,14 +534,21 @@
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_AND_ATT) :: &
b_R_xx_inner_core,b_R_yy_inner_core,b_R_xy_inner_core,b_R_xz_inner_core,b_R_yz_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core,b_epsilondev_xy_inner_core, &
b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
b_eps_trace_over_3_inner_core
+!daniel debug: static
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+! b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core,b_epsilondev_xy_inner_core, &
+! b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core
+!
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+! b_eps_trace_over_3_inner_core
+
! coupling/boundary surfaces
integer :: nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
nspec2D_ymin_inner_core,nspec2D_ymax_inner_core
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90 2013-08-29 07:47:36 UTC (rev 22735)
@@ -0,0 +1,338 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
+! --------------------------------------------------
+!
+! 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
+! August 2013
+!
+! 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.
+!
+!=====================================================================
+
+
+ subroutine update_displacement_Newmark()
+
+! explicit Newmark time scheme with acoustic & elastic domains:
+! (see e.g. Hughes, 1987; Chaljub et al., 2003)
+!
+! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
+! chi_dot(t+delta_t) = chi_dot(t) + 1/2 delta_t chi_dot_dot(t) + 1/2 delta_t chi_dot_dot(t+delta_t)
+! chi_dot_dot(t+delta_t) = 1/M_acoustic( -K_acoustic chi(t+delta) + B_acoustic u(t+delta_t) + f(t+delta_t) )
+!
+! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
+! v(t+delta_t) = v(t) + 1/2 delta_t a(t) + 1/2 delta_t a(t+delta_t)
+! a(t+delta_t) = 1/M_elastic ( -K_elastic u(t+delta) + B_elastic chi_dot_dot(t+delta_t) + f( t+delta_t) )
+!
+! where
+! chi, chi_dot, chi_dot_dot are acoustic (fluid) potentials ( dotted with respect to time)
+! u, v, a are displacement,velocity & acceleration
+! M is mass matrix, K stiffness matrix and B boundary term for acoustic/elastic domains
+! f denotes a source term (acoustic/elastic)
+!
+! note that this stage calculates the predictor terms
+!
+! for
+! potential chi_dot(t+delta) requires + 1/2 delta_t chi_dot_dot(t+delta_t)
+! at a later stage (corrector) once where chi_dot_dot(t+delta) is calculated
+! and similar,
+! velocity v(t+delta_t) requires + 1/2 delta_t a(t+delta_t)
+! at a later stage once where a(t+delta) is calculated
+! also:
+! boundary term B_elastic requires chi_dot_dot(t+delta)
+! thus chi_dot_dot has to be updated first before the elastic boundary term is considered
+
+ use specfem_par
+ use specfem_par_crustmantle
+ use specfem_par_innercore
+ use specfem_par_outercore
+
+ implicit none
+
+ ! updates wavefields
+ if( .not. GPU_MODE) then
+ ! on CPU
+
+ ! Newmark time scheme update
+ ! mantle
+ call update_displ_elastic(NGLOB_CRUST_MANTLE,displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle, &
+ deltat,deltatover2,deltatsqover2)
+ ! outer core
+ call update_displ_acoustic(NGLOB_OUTER_CORE,displ_outer_core,veloc_outer_core,accel_outer_core, &
+ deltat,deltatover2,deltatsqover2)
+ ! inner core
+ call update_displ_elastic(NGLOB_INNER_CORE,displ_inner_core,veloc_inner_core,accel_inner_core, &
+ deltat,deltatover2,deltatsqover2)
+
+ ! backward field
+ if (SIMULATION_TYPE == 3) then
+ ! mantle
+ call update_displ_elastic(NGLOB_CRUST_MANTLE,b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ b_deltat,b_deltatover2,b_deltatsqover2)
+ ! outer core
+ call update_displ_acoustic(NGLOB_OUTER_CORE,b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
+ b_deltat,b_deltatover2,b_deltatsqover2)
+ ! inner core
+ call update_displ_elastic(NGLOB_INNER_CORE,b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
+ b_deltat,b_deltatover2,b_deltatsqover2)
+ endif
+
+ else
+ ! on GPU
+ ! Includes SIM_TYPE 1 & 3
+
+ ! outer core region
+ call it_update_displacement_oc_cuda(Mesh_pointer, &
+ deltat, deltatsqover2, deltatover2, &
+ b_deltat, b_deltatsqover2, b_deltatover2)
+ ! inner core region
+ call it_update_displacement_ic_cuda(Mesh_pointer, &
+ deltat, deltatsqover2, deltatover2, &
+ b_deltat, b_deltatsqover2, b_deltatover2)
+
+ ! crust/mantle region
+ call it_update_displacement_cm_cuda(Mesh_pointer, &
+ deltat, deltatsqover2, deltatover2, &
+ b_deltat, b_deltatsqover2, b_deltatover2)
+ endif
+
+ end subroutine update_displacement_Newmark
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine update_displ_elastic(nglob,displ,veloc,accel, &
+ deltat,deltatover2,deltatsqover2)
+
+ use constants_solver,only: CUSTOM_REAL,NDIM,FORCE_VECTORIZATION_VAL
+
+ implicit none
+
+ integer,intent(in) :: nglob
+ real(kind=CUSTOM_REAL),dimension(NDIM,nglob),intent(inout) :: displ,veloc,accel
+ real(kind=CUSTOM_REAL),intent(in) :: deltat,deltatover2,deltatsqover2
+
+ ! local parameters
+ integer :: i
+
+ ! Newmark time scheme update
+ if(FORCE_VECTORIZATION_VAL) then
+ do i=1,nglob * NDIM
+ displ(i,1) = displ(i,1) + deltat * veloc(i,1) + deltatsqover2 * accel(i,1)
+ veloc(i,1) = veloc(i,1) + deltatover2 * accel(i,1)
+ accel(i,1) = 0._CUSTOM_REAL
+ enddo
+ else
+ do i=1,nglob
+ displ(:,i) = displ(:,i) + deltat * veloc(:,i) + deltatsqover2 * accel(:,i)
+ veloc(:,i) = veloc(:,i) + deltatover2 * accel(:,i)
+ accel(:,i) = 0._CUSTOM_REAL
+ enddo
+ endif
+
+ end subroutine update_displ_elastic
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine update_displ_acoustic(nglob,displ,veloc,accel, &
+ deltat,deltatover2,deltatsqover2)
+
+ use constants,only: CUSTOM_REAL
+
+ implicit none
+
+ integer,intent(in) :: nglob
+ real(kind=CUSTOM_REAL),dimension(nglob),intent(inout) :: displ,veloc,accel
+ real(kind=CUSTOM_REAL),intent(in) :: deltat,deltatover2,deltatsqover2
+
+ ! local parameters
+ integer :: i
+
+ ! Newmark time scheme update
+ do i=1,nglob
+ displ(i) = displ(i) + deltat * veloc(i) + deltatsqover2 * accel(i)
+ veloc(i) = veloc(i) + deltatover2 * accel(i)
+ accel(i) = 0._CUSTOM_REAL
+ enddo
+
+ end subroutine update_displ_acoustic
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine update_veloc_acoustic(NGLOB,veloc_outer_core,accel_outer_core, &
+ deltatover2,rmass_outer_core)
+
+! updates acceleration and velocity in outer core
+
+ use constants_solver,only: CUSTOM_REAL
+
+ implicit none
+
+ integer :: NGLOB
+
+ ! velocity potential
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: veloc_outer_core,accel_outer_core
+
+ ! mass matrix
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmass_outer_core
+
+ real(kind=CUSTOM_REAL) :: deltatover2
+
+ ! local parameters
+ integer :: i
+
+ ! Newmark time scheme
+ ! multiply by the inverse of the mass matrix and update velocity
+
+ do i=1,NGLOB
+ accel_outer_core(i) = accel_outer_core(i)*rmass_outer_core(i)
+ veloc_outer_core(i) = veloc_outer_core(i) + deltatover2*accel_outer_core(i)
+ enddo
+
+ end subroutine update_veloc_acoustic
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine update_accel_elastic(NGLOB,NGLOB_XY,veloc_crust_mantle,accel_crust_mantle, &
+ two_omega_earth, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NCHUNKS_VAL,ABSORBING_CONDITIONS)
+
+! updates acceleration in crust/mantle region
+
+ use constants_solver,only: CUSTOM_REAL,NDIM
+
+ implicit none
+
+ integer :: NGLOB,NGLOB_XY,NCHUNKS_VAL
+
+ ! velocity & acceleration
+ ! crust/mantle region
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc_crust_mantle,accel_crust_mantle
+
+ ! mass matrices
+ !
+ ! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmassz_crust_mantle
+
+ real(kind=CUSTOM_REAL) :: two_omega_earth
+
+ logical :: ABSORBING_CONDITIONS
+
+ ! local parameters
+ integer :: i
+
+ ! updates acceleration w/ rotation in crust/mantle region only
+
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+
+ do i=1,NGLOB
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
+ - 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
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+
+ endif
+
+ end subroutine update_accel_elastic
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine update_veloc_elastic(NGLOB_CM,veloc_crust_mantle,accel_crust_mantle, &
+ NGLOB_IC,veloc_inner_core,accel_inner_core, &
+ deltatover2,two_omega_earth,rmass_inner_core)
+
+! updates velocity in crust/mantle region, and acceleration and velocity in inner core
+
+ use constants_solver,only: CUSTOM_REAL,NDIM
+
+ implicit none
+
+ integer :: NGLOB_CM,NGLOB_IC
+
+ ! acceleration & velocity
+ ! crust/mantle region
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CM) :: veloc_crust_mantle,accel_crust_mantle
+ ! inner core
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_IC) :: veloc_inner_core,accel_inner_core
+
+ ! mass matrix
+ real(kind=CUSTOM_REAL), dimension(NGLOB_IC) :: rmass_inner_core
+
+ real(kind=CUSTOM_REAL) :: deltatover2,two_omega_earth
+
+ ! local parameters
+ integer :: i
+
+ ! Newmark time scheme:
+ !
+ ! note:
+ ! - crust/mantle region
+ ! needs only velocity corrector terms
+ ! (acceleration already updated before)
+ ! - inner core region
+ ! needs both, acceleration update & velocity corrector terms
+
+ ! mantle
+ do i=1,NGLOB_CM
+ veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
+ enddo
+
+ ! inner core
+ do i=1,NGLOB_IC
+ 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)
+
+ veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
+ enddo
+
+ end subroutine update_veloc_elastic
More information about the CIG-COMMITS
mailing list