[cig-commits] r20103 - in seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER: in_data_files src/cuda src/decompose_mesh_SCOTCH src/generate_databases src/meshfem3D src/shared src/specfem3D utils
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sun May 13 13:21:07 PDT 2012
Author: danielpeter
Date: 2012-05-13 13:21:06 -0700 (Sun, 13 May 2012)
New Revision: 20103
Added:
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/update_headers_change_word_f90.pl
Modified:
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/in_data_files/Par_file.in
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_model.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/memory_eval.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_cascadia.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_prem.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_socal.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_default.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_external_values.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_gll.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_salton_trough.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_tomography.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/check_mesh_quality.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/compute_parameters.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_visual_files.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/store_boundaries.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/combine_vol_data.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_topo_bathy_file.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/smooth_vol_data.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/assemble_MPI_vector.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_po.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_poroelastic_el.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_fluid.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/create_color_image.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/noise_tomography.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/write_movie_output.f90
Log:
cleans up, includes current trunk routines
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/in_data_files/Par_file.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/in_data_files/Par_file.in 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/in_data_files/Par_file.in 2012-05-13 20:21:06 UTC (rev 20103)
@@ -15,6 +15,15 @@
NSTEP = 1000
DT = 0.05d0
+# models:
+# available options are:
+# default (model parameters described by mesh properties)
+# 1D models available are:
+# 1d_prem,1d_socal,1d_cascadia
+# 3D models available are:
+# aniso,external,gll,salton_trough,tomo
+MODEL = default
+
# parameters describing the model
OCEANS = .false.
TOPOGRAPHY = .false.
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-05-13 20:21:06 UTC (rev 20103)
@@ -450,6 +450,7 @@
// Gets rank number of MPI process
int myrank = *myrank_f;
+/*
// cuda initialization (needs -lcuda library)
// note: cuInit initializes the driver API.
// it is needed for any following CUDA driver API function call (format cuFUNCTION(..) )
@@ -477,12 +478,14 @@
fprintf(stderr,"Compute capability should be at least 1.3, got: %d.%d \nexiting...\n",major,minor);
exit_on_error("CUDA Compute capability major number should be at least 1.3\n");
}
+*/
// note: from here on we use the runtime API ...
+
// Gets number of GPU devices
int device_count = 0;
cudaGetDeviceCount(&device_count);
- exit_on_cuda_error("CUDA runtime cudaGetDeviceCount: check if driver and runtime libraries work together\nexiting...\n");
+ exit_on_cuda_error("CUDA runtime error: cudaGetDeviceCount failed\ncheck if driver and runtime libraries work together\nexiting...\n");
// returns device count to fortran
if (device_count == 0) exit_on_error("CUDA runtime error: there is no device supporting CUDA\n");
@@ -548,16 +551,6 @@
fprintf(fp,"deviceOverlap: FALSE\n");
}
- // make sure that the device has compute capability >= 1.3
- //if (deviceProp.major < 1){
- // fprintf(stderr,"Compute capability major number should be at least 1, exiting...\n\n");
- // exit_on_error("CUDA Compute capability major number should be at least 1");
- //}
- //if (deviceProp.major == 1 && deviceProp.minor < 3){
- // fprintf(stderr,"Compute capability should be at least 1.3, exiting...\n");
- // exit_on_error("CUDA Compute capability major number should be at least 1.3");
- //}
-
// outputs initial memory infos via cudaMemGetInfo()
double free_db,used_db,total_db;
get_free_memory(&free_db,&used_db,&total_db);
@@ -566,6 +559,17 @@
fclose(fp);
}
+
+ // make sure that the device has compute capability >= 1.3
+ if (deviceProp.major < 1){
+ fprintf(stderr,"Compute capability major number should be at least 1, exiting...\n\n");
+ exit_on_error("CUDA Compute capability major number should be at least 1\n");
+ }
+ if (deviceProp.major == 1 && deviceProp.minor < 3){
+ fprintf(stderr,"Compute capability should be at least 1.3, exiting...\n");
+ exit_on_error("CUDA Compute capability major number should be at least 1.3\n");
+ }
+
}
/* ----------------------------------------------------------------------------------------------- */
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -41,7 +41,7 @@
integer, dimension(:,:), allocatable :: elmnts
integer, dimension(:,:), allocatable :: mat
integer, dimension(:), allocatable :: part
-
+
integer :: nnodes
double precision, dimension(:,:), allocatable :: nodes_coords
@@ -121,7 +121,7 @@
character(len=256) :: line
logical :: use_poroelastic_file
integer(long) :: nspec_long
-
+
! sets number of nodes per element
ngnod = esize
@@ -158,10 +158,10 @@
print*,'bit size fortran: ',bit_size(nspec)
stop 'error number of elements too large'
endif
-
+
! sets number of elements (integer 4-byte)
nspec = nspec_long
-
+
allocate(elmnts(esize,nspec),stat=ier)
if( ier /= 0 ) stop 'error allocating array elmnts'
do ispec = 1, nspec
@@ -810,7 +810,7 @@
if (ier /= 0) then
stop 'ERROR : MAIN : Cannot destroy strat'
endif
-
+
! re-partitioning puts poroelastic-elastic coupled elements into same partition
! integer :: nfaces_coupled
! integer, dimension(:,:), pointer :: faces_coupled
@@ -826,7 +826,7 @@
nspec2D_moho,ibelm_moho,nodes_ibelm_moho )
- ! local number of each element for each partition
+ ! local number of each element for each partition
call build_glob2loc_elmnts(nspec, part, glob2loc_elmnts,nparts)
! local number of each node for each partition
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in 2012-05-13 20:21:06 UTC (rev 20103)
@@ -200,12 +200,7 @@
$O/parallel.o: ${SHARED}/constants.h ${SHARED}/parallel.f90
${MPIFCCOMPILE_CHECK} -c -o $O/parallel.o ${SHARED}/parallel.f90
-$O/model_external_values.o: ${SHARED}/constants.h model_external_values.f90
- ${FCCOMPILE_CHECK} -c -o $O/model_external_values.o model_external_values.f90
-$O/model_tomography.o: ${SHARED}/constants.h model_tomography.f90
- ${FCCOMPILE_CHECK} -c -o $O/model_tomography.o model_tomography.f90
-
###
### serial compilation without optimization
###
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_mass_matrices.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_mass_matrices.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -160,7 +160,7 @@
integer :: ier
real(kind=CUSTOM_REAL) :: xloc,yloc,loc_elevation
-
+
! creates ocean load mass matrix
if(OCEANS) then
@@ -191,17 +191,17 @@
! compute local height of oceans
if( TOPOGRAPHY ) then
-
+
! takes elevation from topography file
xloc = xstore_dummy(iglobnum)
yloc = ystore_dummy(iglobnum)
-
+
call get_topo_bathy_elevation(xloc,yloc,loc_elevation, &
itopo_bathy,NX_TOPO,NY_TOPO, &
UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION)
-
+
elevation = dble(loc_elevation)
-
+
else
! takes elevation from z-coordinate of mesh point
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -33,7 +33,7 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: ystore_dummy
real(kind=CUSTOM_REAL), dimension(:), allocatable :: zstore_dummy
integer :: nglob_dummy
-
+
! Gauss-Lobatto-Legendre points and weights of integration
double precision, dimension(:), allocatable :: xigll,yigll,zigll,wxgll,wygll,wzgll
@@ -148,7 +148,12 @@
logical :: ACOUSTIC_SIMULATION,ELASTIC_SIMULATION,POROELASTIC_SIMULATION
+ ! mesh coloring
+ integer :: num_colors_outer_acoustic,num_colors_inner_acoustic
+ integer, dimension(:), allocatable :: num_elem_colors_acoustic
+ integer :: num_colors_outer_elastic,num_colors_inner_elastic
+ integer, dimension(:), allocatable :: num_elem_colors_elastic
end module create_regions_mesh_ext_par
!
@@ -182,15 +187,10 @@
ATTENUATION,USE_OLSEN_ATTENUATION, &
UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
NX_TOPO,NY_TOPO,itopo_bathy, &
- nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho
+ nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho
use create_regions_mesh_ext_par
implicit none
-! moho (optional)
- integer :: nspec2D_moho_ext
- integer, dimension(nspec2D_moho_ext) :: ibelm_moho
- integer, dimension(4,nspec2D_moho_ext) :: nodes_ibelm_moho
-
! local parameters
! static memory size needed by the solver
double precision :: static_memory_size
@@ -340,6 +340,13 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
ibool,SAVE_MESH_FILES)
+! colors mesh if requested
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...element mesh coloring '
+ endif
+ call crm_setup_color_perm(myrank,nspec,nglob,ibool,ANISOTROPY,SAVE_MESH_FILES)
+
! saves the binary mesh files
call sync_all()
if( myrank == 0) then
@@ -384,7 +391,11 @@
nspec_outer_acoustic,nspec_outer_elastic,nspec_outer_poroelastic, &
num_phase_ispec_acoustic,phase_ispec_inner_acoustic, &
num_phase_ispec_elastic,phase_ispec_inner_elastic, &
- num_phase_ispec_poroelastic,phase_ispec_inner_poroelastic)
+ num_phase_ispec_poroelastic,phase_ispec_inner_poroelastic, &
+ num_colors_outer_acoustic,num_colors_inner_acoustic, &
+ num_elem_colors_acoustic, &
+ num_colors_outer_elastic,num_colors_inner_elastic, &
+ num_elem_colors_elastic)
! saves moho surface
if( SAVE_MOHO_MESH ) then
@@ -1359,3 +1370,649 @@
end subroutine crm_setup_inner_outer_elemnts
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine crm_setup_color_perm(myrank,nspec,nglob,ibool,ANISOTROPY,SAVE_MESH_FILES)
+
+! sets up mesh coloring and permutes elements
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ integer :: myrank,nspec,nglob
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+ logical :: ANISOTROPY,SAVE_MESH_FILES
+
+ ! local parameters
+ integer, dimension(:), allocatable :: perm
+ integer :: ier
+
+ ! user output
+ if(myrank == 0) then
+ write(IMAIN,*) ' use coloring = ',USE_MESH_COLORING_GPU
+ endif
+
+ ! initializes
+ num_colors_outer_acoustic = 0
+ num_colors_inner_acoustic = 0
+ num_colors_outer_elastic = 0
+ num_colors_inner_elastic = 0
+
+ ! mesh coloring
+ if( USE_MESH_COLORING_GPU ) then
+
+ ! creates coloring of elements
+ allocate(perm(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating temporary perm array'
+ perm(:) = 0
+
+ ! acoustic domains
+ if( ACOUSTIC_SIMULATION ) then
+ if( myrank == 0) then
+ write(IMAIN,*) ' acoustic domains: '
+ endif
+ call crm_setup_color(myrank,nspec,nglob,ibool,perm, &
+ ispec_is_acoustic,1, &
+ num_phase_ispec_acoustic,phase_ispec_inner_acoustic, &
+ SAVE_MESH_FILES)
+ else
+ ! allocates dummy arrays
+ allocate(num_elem_colors_acoustic(num_colors_outer_acoustic + num_colors_inner_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_acoustic array'
+ endif
+
+ ! elastic domains
+ if( ELASTIC_SIMULATION ) then
+ if( myrank == 0) then
+ write(IMAIN,*) ' elastic domains: '
+ endif
+ call crm_setup_color(myrank,nspec,nglob,ibool,perm, &
+ ispec_is_elastic,2, &
+ num_phase_ispec_elastic,phase_ispec_inner_elastic, &
+ SAVE_MESH_FILES)
+ else
+ allocate(num_elem_colors_elastic(num_colors_outer_elastic + num_colors_inner_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_elastic array'
+ endif
+
+ ! checks: after all domains are done
+ if(minval(perm) /= 1) &
+ call exit_MPI(myrank, 'minval(perm) should be 1')
+ if(maxval(perm) /= max(num_phase_ispec_acoustic,num_phase_ispec_elastic)) &
+ call exit_MPI(myrank, 'maxval(perm) should be max(num_phase_ispec_..)')
+
+ ! sorts array according to permutation
+ call sync_all()
+ if(myrank == 0) then
+ write(IMAIN,*) ' mesh permutation:'
+ endif
+ call crm_setup_permutation(myrank,nspec,nglob,ibool,ANISOTROPY,perm, &
+ SAVE_MESH_FILES)
+
+ deallocate(perm)
+
+ else
+
+ ! allocates dummy arrays
+ allocate(num_elem_colors_acoustic(num_colors_outer_acoustic + num_colors_inner_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_acoustic array'
+ allocate(num_elem_colors_elastic(num_colors_outer_elastic + num_colors_inner_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_elastic array'
+
+ endif ! USE_MESH_COLORING_GPU
+
+ end subroutine crm_setup_color_perm
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine crm_setup_color(myrank,nspec,nglob,ibool,perm, &
+ ispec_is_d,idomain, &
+ num_phase_ispec_d,phase_ispec_inner_d, &
+ SAVE_MESH_FILES)
+
+! sets up mesh coloring
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ integer :: myrank,nspec,nglob
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+ integer, dimension(nspec) :: perm
+
+ ! wrapper array for ispec is in domain:
+ ! idomain acoustic == 1 or elastic == 2
+ integer :: idomain
+ logical, dimension(nspec) :: ispec_is_d
+ integer :: num_phase_ispec_d
+ integer, dimension(num_phase_ispec_d,2) :: phase_ispec_inner_d
+
+ logical :: SAVE_MESH_FILES
+
+ ! local parameters
+ ! added for color permutation
+ integer :: nb_colors_outer_elements,nb_colors_inner_elements
+ integer, dimension(:), allocatable :: num_of_elems_in_this_color
+ integer, dimension(:), allocatable :: color
+ integer, dimension(:), allocatable :: first_elem_number_in_this_color
+ logical, dimension(:), allocatable :: is_on_a_slice_edge
+
+ integer :: nspec_outer,nspec_inner,nspec_domain
+ integer :: nspec_outer_min_global,nspec_outer_max_global
+ integer :: nb_colors,nb_colors_min,nb_colors_max
+
+ integer :: icolor,ispec,ispec_counter
+ integer :: ispec_inner,ispec_outer
+ integer :: ier
+
+ character(len=2),dimension(2) :: str_domain = (/ "ac", "el" /)
+ character(len=256) :: filename
+
+ !!!! David Michea: detection of the edges, coloring and permutation separately
+
+ ! implement mesh coloring for GPUs if needed, to create subsets of disconnected elements
+ ! to remove dependencies and the need for atomic operations in the sum of
+ ! elemental contributions in the solver
+
+ ! allocates temporary array with colors
+ allocate(color(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating temporary color array'
+ allocate(first_elem_number_in_this_color(MAX_NUMBER_OF_COLORS + 1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating first_elem_number_in_this_color array'
+
+ ! flags for elements on outer rims
+ ! opposite to what is stored in ispec_is_inner
+ allocate(is_on_a_slice_edge(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating is_on_a_slice_edge array'
+ do ispec = 1,nspec
+ is_on_a_slice_edge(ispec) = .not. ispec_is_inner(ispec)
+ enddo
+
+ ! fast element coloring scheme
+ call get_perm_color_faster(is_on_a_slice_edge,ispec_is_d, &
+ ibool,perm,color, &
+ nspec,nglob, &
+ nb_colors_outer_elements,nb_colors_inner_elements, &
+ nspec_outer,nspec_inner,nspec_domain, &
+ first_elem_number_in_this_color, &
+ myrank)
+
+ ! for the last color, the next color is fictitious and its first (fictitious) element number is nspec + 1
+ first_elem_number_in_this_color(nb_colors_outer_elements + nb_colors_inner_elements + 1) &
+ = nspec_domain + 1
+
+ allocate(num_of_elems_in_this_color(nb_colors_outer_elements + nb_colors_inner_elements),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_of_elems_in_this_color array'
+
+ num_of_elems_in_this_color(:) = 0
+ do icolor = 1, nb_colors_outer_elements + nb_colors_inner_elements
+ num_of_elems_in_this_color(icolor) = first_elem_number_in_this_color(icolor+1) - first_elem_number_in_this_color(icolor)
+ enddo
+
+ ! check that the sum of all the numbers of elements found in each color is equal
+ ! to the total number of elements in the mesh
+ if(sum(num_of_elems_in_this_color) /= nspec_domain) then
+ print *,'error number of elements in this color:',idomain
+ print *,'rank: ',myrank,' nspec = ',nspec_domain
+ print *,' total number of elements in all the colors of the mesh = ', &
+ sum(num_of_elems_in_this_color)
+ call exit_MPI(myrank, 'incorrect total number of elements in all the colors of the mesh')
+ endif
+
+ ! check that the sum of all the numbers of elements found in each color for the outer elements is equal
+ ! to the total number of outer elements found in the mesh
+ if(sum(num_of_elems_in_this_color(1:nb_colors_outer_elements)) /= nspec_outer) then
+ print *,'error number of outer elements in this color:',idomain
+ print *,'rank: ',myrank,' nspec_outer = ',nspec_outer
+ print*,'nb_colors_outer_elements = ',nb_colors_outer_elements
+ print *,'total number of elements in all the colors of the mesh for outer elements = ', &
+ sum(num_of_elems_in_this_color(1:nb_colors_outer_elements))
+ call exit_MPI(myrank, 'incorrect total number of elements in all the colors of the mesh for outer elements')
+ endif
+
+ ! debug: file output
+ if( SAVE_MESH_FILES ) then
+ filename = prname(1:len_trim(prname))//'color_'//str_domain(idomain)
+ call write_VTK_data_elem_i(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ color,filename)
+ endif
+
+ deallocate(first_elem_number_in_this_color)
+ deallocate(is_on_a_slice_edge)
+ deallocate(color)
+
+ ! debug: no mesh coloring, only creates dummy coloring arrays
+ if( .false. ) then
+ nb_colors_outer_elements = 0
+ nb_colors_inner_elements = 0
+ ispec_counter = 0
+
+ ! first generate all the outer elements
+ do ispec = 1,nspec
+ if( ispec_is_d(ispec) ) then
+ if( ispec_is_inner(ispec) .eqv. .false. ) then
+ ispec_counter = ispec_counter + 1
+ perm(ispec) = ispec_counter
+ endif
+ endif
+ enddo
+
+ ! store total number of outer elements
+ nspec_outer = ispec_counter
+
+ ! only single color
+ if(nspec_outer > 0 ) nb_colors_outer_elements = 1
+
+ ! then generate all the inner elements
+ do ispec = 1,nspec
+ if( ispec_is_d(ispec) ) then
+ if( ispec_is_inner(ispec) .eqv. .true. ) then
+ ispec_counter = ispec_counter + 1
+ perm(ispec) = ispec_counter - nspec_outer ! starts again at 1
+ endif
+ endif
+ enddo
+ nspec_inner = ispec_counter - nspec_outer
+
+ ! only single color
+ if(nspec_inner > 0 ) nb_colors_inner_elements = 1
+
+ allocate(num_of_elems_in_this_color(nb_colors_outer_elements + nb_colors_inner_elements),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_of_elems_in_this_color array'
+
+ if( nspec_outer > 0 ) num_of_elems_in_this_color(1) = nspec_outer
+ if( nspec_inner > 0 ) num_of_elems_in_this_color(2) = nspec_inner
+ endif ! debug
+
+ ! debug: saves mesh coloring numbers into files
+ if( SAVE_MESH_FILES ) then
+ ! debug file output
+ filename = prname(1:len_trim(prname))//'num_of_elems_in_this_color_'//str_domain(idomain)//'.dat'
+ open(unit=99,file=trim(filename),status='unknown',iostat=ier)
+ if( ier /= 0 ) stop 'error opening num_of_elems_in_this_color file'
+ ! number of colors for outer elements
+ write(99,*) nb_colors_outer_elements
+ ! number of colors for inner elements
+ write(99,*) nb_colors_inner_elements
+ ! number of elements in each color
+ ! outer elements
+ do icolor = 1, nb_colors_outer_elements + nb_colors_inner_elements
+ write(99,*) num_of_elems_in_this_color(icolor)
+ enddo
+ close(99)
+ endif
+
+ ! sets up domain coloring arrays
+ select case(idomain)
+ case( 1 )
+ ! acoustic domains
+ num_colors_outer_acoustic = nb_colors_outer_elements
+ num_colors_inner_acoustic = nb_colors_inner_elements
+
+ allocate(num_elem_colors_acoustic(num_colors_outer_acoustic + num_colors_inner_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_acoustic array'
+
+ num_elem_colors_acoustic(:) = num_of_elems_in_this_color(:)
+
+ case( 2 )
+ ! elastic domains
+ num_colors_outer_elastic = nb_colors_outer_elements
+ num_colors_inner_elastic = nb_colors_inner_elements
+
+ allocate(num_elem_colors_elastic(num_colors_outer_elastic + num_colors_inner_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_elastic array'
+
+ num_elem_colors_elastic(:) = num_of_elems_in_this_color(:)
+
+ case default
+ stop 'error idomain not recognized'
+ end select
+
+ ! sets up elements for loops in simulations
+ ispec_inner = 0
+ ispec_outer = 0
+ do ispec = 1, nspec
+ ! only elements in this domain
+ if( ispec_is_d(ispec) ) then
+
+ ! sets phase_ispec arrays with ordering of elements
+ if( ispec_is_inner(ispec) .eqv. .false. ) then
+ ! outer elements
+ ispec_outer = perm(ispec)
+
+ ! checks
+ if( ispec_outer < 1 .or. ispec_outer > num_phase_ispec_d ) then
+ print*,'error outer permutation:',idomain
+ print*,'rank:',myrank,' ispec_inner = ',ispec_outer
+ print*,'num_phase_ispec_d = ',num_phase_ispec_d
+ call exit_MPI(myrank,'error outer acoustic permutation')
+ endif
+
+ phase_ispec_inner_d(ispec_outer,1) = ispec
+
+ else
+ ! inner elements
+ ispec_inner = perm(ispec)
+
+ ! checks
+ if( ispec_inner < 1 .or. ispec_inner > num_phase_ispec_d ) then
+ print*,'error inner permutation:',idomain
+ print*,'rank:',myrank,' ispec_inner = ',ispec_inner
+ print*,'num_phase_ispec_d = ',num_phase_ispec_d
+ call exit_MPI(myrank,'error inner acoustic permutation')
+ endif
+
+ phase_ispec_inner_d(ispec_inner,2) = ispec
+
+ endif
+ endif
+ enddo
+
+ ! total number of colors
+ nb_colors = nb_colors_inner_elements + nb_colors_outer_elements
+ call min_all_i(nb_colors,nb_colors_min)
+ call max_all_i(nb_colors,nb_colors_max)
+
+ ! user output
+ call min_all_i(nspec_outer,nspec_outer_min_global)
+ call max_all_i(nspec_outer,nspec_outer_max_global)
+ call min_all_i(nspec_outer,nspec_outer_min_global)
+ call max_all_i(nspec_outer,nspec_outer_max_global)
+ if(myrank == 0) then
+ write(IMAIN,*) ' colors min = ',nb_colors_min
+ write(IMAIN,*) ' colors max = ',nb_colors_max
+ write(IMAIN,*) ' outer elements: min = ',nspec_outer_min_global
+ write(IMAIN,*) ' outer elements: max = ',nspec_outer_max_global
+ endif
+
+ ! debug: outputs permutation array as vtk file
+ !if( SAVE_MESH_FILES ) then
+ ! filename = prname(1:len_trim(prname))//'perm_'//str_domain(idomain)
+ ! call write_VTK_data_elem_i(nspec,nglob, &
+ ! xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ ! perm,filename)
+ !endif
+
+ deallocate(num_of_elems_in_this_color)
+
+ end subroutine crm_setup_color
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine crm_setup_permutation(myrank,nspec,nglob,ibool,ANISOTROPY,perm, &
+ SAVE_MESH_FILES)
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ integer :: myrank,nspec,nglob
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+ logical :: ANISOTROPY
+
+ integer, dimension(nspec),intent(inout) :: perm
+
+ logical :: SAVE_MESH_FILES
+
+ ! local parameters
+ ! added for sorting
+ integer, dimension(:,:,:,:), allocatable :: temp_array_int
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: temp_array_real
+ logical, dimension(:), allocatable :: temp_array_logical_1D
+
+ integer, dimension(:), allocatable :: temp_perm_global
+ logical, dimension(:), allocatable :: mask_global
+
+ integer :: icolor,icounter,ispec,ielem,ier,i
+ integer :: iface,old_ispec,new_ispec
+
+ character(len=256) :: filename
+
+ ! sorts array according to permutation
+ allocate(temp_perm_global(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error temp_perm_global array'
+
+ ! global ordering
+ temp_perm_global(:) = 0
+ icounter = 0
+
+ ! fills global permutation array
+ ! starts with elastic elements
+ if( ELASTIC_SIMULATION ) then
+ ! first outer elements coloring
+ ! phase element counter
+ ielem = 0
+ do icolor = 1,num_colors_outer_elastic
+ ! loops through elements
+ do i = 1,num_elem_colors_elastic(icolor)
+ ielem = ielem + 1
+ ispec = phase_ispec_inner_elastic(ielem,1) ! 1 <-- first phase, outer elements
+ ! reorders elements
+ icounter = icounter + 1
+ temp_perm_global(ispec) = icounter
+ ! resets to new order
+ phase_ispec_inner_elastic(ielem,1) = icounter
+ enddo
+ enddo
+ ! inner elements coloring
+ ielem = 0
+ do icolor = num_colors_outer_elastic+1,num_colors_outer_elastic+num_colors_inner_elastic
+ ! loops through elements
+ do i = 1,num_elem_colors_elastic(icolor)
+ ielem = ielem + 1
+ ispec = phase_ispec_inner_elastic(ielem,2) ! 2 <-- second phase, inner elements
+ ! reorders elements
+ icounter = icounter + 1
+ temp_perm_global(ispec) = icounter
+ ! resets to new order
+ phase_ispec_inner_elastic(ielem,2) = icounter
+ enddo
+ enddo
+ endif
+
+ ! continues with acoustic elements
+ if( ACOUSTIC_SIMULATION ) then
+ ! first outer elements coloring
+ ! phase element counter
+ ielem = 0
+ do icolor = 1,num_colors_outer_acoustic
+ ! loops through elements
+ do i = 1,num_elem_colors_acoustic(icolor)
+ ielem = ielem + 1
+ ispec = phase_ispec_inner_acoustic(ielem,1) ! 1 <-- first phase, outer elements
+ ! reorders elements
+ icounter = icounter + 1
+ temp_perm_global(ispec) = icounter
+ ! resets to new order
+ phase_ispec_inner_acoustic(ielem,1) = icounter
+ enddo
+ enddo
+ ! inner elements coloring
+ ielem = 0
+ do icolor = num_colors_outer_acoustic+1,num_colors_outer_acoustic+num_colors_inner_acoustic
+ ! loops through elements
+ do i = 1,num_elem_colors_acoustic(icolor)
+ ielem = ielem + 1
+ ispec = phase_ispec_inner_acoustic(ielem,2) ! 2 <-- second phase, inner elements
+ ! reorders elements
+ icounter = icounter + 1
+ temp_perm_global(ispec) = icounter
+ ! resets to new order
+ phase_ispec_inner_acoustic(ielem,2) = icounter
+ enddo
+ enddo
+ endif
+
+ ! checks
+ if( icounter /= nspec ) then
+ print*,'error temp perm: ',icounter,nspec
+ stop 'error temporary global permutation incomplete'
+ endif
+
+ ! checks perm entries
+ if(minval(temp_perm_global) /= 1) call exit_MPI(myrank, 'minval(temp_perm_global) should be 1')
+ if(maxval(temp_perm_global) /= nspec) call exit_MPI(myrank, 'maxval(temp_perm_global) should be nspec')
+
+ ! checks if every element was uniquely set
+ allocate(mask_global(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating temporary mask_global'
+ mask_global(:) = .false.
+ icounter = 0 ! counts permutations
+ do ispec = 1, nspec
+ new_ispec = temp_perm_global(ispec)
+ ! checks bounds
+ if( new_ispec < 1 .or. new_ispec > nspec ) call exit_MPI(myrank,'error temp_perm_global ispec bounds')
+ ! checks if already set
+ if( mask_global(new_ispec) ) then
+ print*,'error temp_perm_global:',ispec,new_ispec,'element already set'
+ call exit_MPI(myrank,'error global permutation')
+ else
+ mask_global(new_ispec) = .true.
+ endif
+ ! counts permutations
+ if( new_ispec /= ispec ) icounter = icounter + 1
+ enddo
+
+ ! checks number of set elements
+ if( count(mask_global(:)) /= nspec ) then
+ print*,'error temp_perm_global:',count(mask_global(:)),nspec,'permutation incomplete'
+ call exit_MPI(myrank,'error global permutation incomplete')
+ endif
+ deallocate(mask_global)
+
+ ! user output
+ if(myrank == 0) then
+ write(IMAIN,*) ' number of permutations = ',icounter
+ endif
+
+ ! outputs permutation array as vtk file
+ if( SAVE_MESH_FILES ) then
+ filename = prname(1:len_trim(prname))//'perm_global'
+ call write_VTK_data_elem_i(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ temp_perm_global,filename)
+ endif
+
+ ! store as new permutation
+ perm(:) = temp_perm_global(:)
+ deallocate(temp_perm_global)
+
+ ! permutes all required mesh arrays according to new ordering
+
+ ! permutation of ibool
+ allocate(temp_array_int(NGLLX,NGLLY,NGLLZ,nspec))
+ call permute_elements_integer(ibool,temp_array_int,perm,nspec)
+ deallocate(temp_array_int)
+
+ ! element domain flags
+ allocate(temp_array_logical_1D(nspec))
+ call permute_elements_logical1D(ispec_is_acoustic,temp_array_logical_1D,perm,nspec)
+ call permute_elements_logical1D(ispec_is_elastic,temp_array_logical_1D,perm,nspec)
+ call permute_elements_logical1D(ispec_is_poroelastic,temp_array_logical_1D,perm,nspec)
+ call permute_elements_logical1D(ispec_is_inner,temp_array_logical_1D,perm,nspec)
+ deallocate(temp_array_logical_1D)
+
+ ! mesh arrays
+ allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec))
+ call permute_elements_real(xixstore,temp_array_real,perm,nspec)
+ call permute_elements_real(xiystore,temp_array_real,perm,nspec)
+ call permute_elements_real(xizstore,temp_array_real,perm,nspec)
+ call permute_elements_real(etaxstore,temp_array_real,perm,nspec)
+ call permute_elements_real(etaystore,temp_array_real,perm,nspec)
+ call permute_elements_real(etazstore,temp_array_real,perm,nspec)
+ call permute_elements_real(gammaxstore,temp_array_real,perm,nspec)
+ call permute_elements_real(gammaystore,temp_array_real,perm,nspec)
+ call permute_elements_real(gammazstore,temp_array_real,perm,nspec)
+ call permute_elements_real(jacobianstore,temp_array_real,perm,nspec)
+
+ call permute_elements_real(kappastore,temp_array_real,perm,nspec)
+ call permute_elements_real(mustore,temp_array_real,perm,nspec)
+
+ ! acoustic arrays
+ if( ACOUSTIC_SIMULATION ) then
+ call permute_elements_real(rhostore,temp_array_real,perm,nspec)
+ endif
+
+ ! elastic arrays
+ if( ELASTIC_SIMULATION ) then
+ call permute_elements_real(rho_vp,temp_array_real,perm,nspec)
+ call permute_elements_real(rho_vs,temp_array_real,perm,nspec)
+ if( ANISOTROPY ) then
+ call permute_elements_real(c11store,temp_array_real,perm,nspec)
+ call permute_elements_real(c12store,temp_array_real,perm,nspec)
+ call permute_elements_real(c13store,temp_array_real,perm,nspec)
+ call permute_elements_real(c14store,temp_array_real,perm,nspec)
+ call permute_elements_real(c15store,temp_array_real,perm,nspec)
+ call permute_elements_real(c16store,temp_array_real,perm,nspec)
+ call permute_elements_real(c22store,temp_array_real,perm,nspec)
+ call permute_elements_real(c23store,temp_array_real,perm,nspec)
+ call permute_elements_real(c24store,temp_array_real,perm,nspec)
+ call permute_elements_real(c25store,temp_array_real,perm,nspec)
+ call permute_elements_real(c33store,temp_array_real,perm,nspec)
+ call permute_elements_real(c34store,temp_array_real,perm,nspec)
+ call permute_elements_real(c35store,temp_array_real,perm,nspec)
+ call permute_elements_real(c36store,temp_array_real,perm,nspec)
+ call permute_elements_real(c44store,temp_array_real,perm,nspec)
+ call permute_elements_real(c45store,temp_array_real,perm,nspec)
+ call permute_elements_real(c46store,temp_array_real,perm,nspec)
+ call permute_elements_real(c55store,temp_array_real,perm,nspec)
+ call permute_elements_real(c56store,temp_array_real,perm,nspec)
+ call permute_elements_real(c66store,temp_array_real,perm,nspec)
+ endif
+ endif
+ deallocate(temp_array_real)
+
+ ! boundary surface
+ if( num_abs_boundary_faces > 0 ) then
+ do iface = 1,num_abs_boundary_faces
+ old_ispec = abs_boundary_ispec(iface)
+ new_ispec = perm(old_ispec)
+ abs_boundary_ispec(iface) = new_ispec
+ enddo
+ endif
+
+ ! free surface
+ if( num_free_surface_faces > 0 ) then
+ do iface = 1,num_free_surface_faces
+ old_ispec = free_surface_ispec(iface)
+ new_ispec = perm(old_ispec)
+ free_surface_ispec(iface) = new_ispec
+ enddo
+ endif
+
+ ! coupling surface
+ if( num_coupling_ac_el_faces > 0 ) then
+ do iface = 1,num_coupling_ac_el_faces
+ old_ispec = coupling_ac_el_ispec(iface)
+ new_ispec = perm(old_ispec)
+ coupling_ac_el_ispec(iface) = new_ispec
+ enddo
+ endif
+
+ ! moho surface
+ if( NSPEC2D_MOHO > 0 ) then
+ allocate(temp_array_logical_1D(nspec))
+ call permute_elements_logical1D(is_moho_top,temp_array_logical_1D,perm,nspec)
+ call permute_elements_logical1D(is_moho_bot,temp_array_logical_1D,perm,nspec)
+ deallocate(temp_array_logical_1D)
+ do iface = 1,NSPEC2D_MOHO
+ ! top
+ old_ispec = ibelm_moho_top(iface)
+ new_ispec = perm(old_ispec)
+ ibelm_moho_top(iface) = new_ispec
+ ! bottom
+ old_ispec = ibelm_moho_bot(iface)
+ new_ispec = perm(old_ispec)
+ ibelm_moho_bot(iface) = new_ispec
+ enddo
+ endif
+
+ end subroutine crm_setup_permutation
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/generate_databases.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/generate_databases.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -295,7 +295,7 @@
! flag for noise simulation
integer :: NOISE_TOMOGRAPHY
integer :: IMODEL
-
+
end module generate_databases_par
!
@@ -418,7 +418,7 @@
write(IMAIN,*) 'Shape functions defined by NGNOD = ',NGNOD,' control nodes'
write(IMAIN,*) 'Surface shape functions defined by NGNOD2D = ',NGNOD2D,' control nodes'
write(IMAIN,*)
-
+
write(IMAIN,'(a)',advance='no') ' velocity model: '
select case(IMODEL)
case( IMODEL_DEFAULT )
@@ -438,7 +438,7 @@
case( IMODEL_USER_EXTERNAL )
write(IMAIN,'(a)',advance='yes') ' external'
end select
-
+
write(IMAIN,*)
endif
@@ -588,8 +588,8 @@
open(unit=IIN,file=prname(1:len_trim(prname))//'Database', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
- write(IMAIN,*) 'error opening file: ',prname(1:len_trim(prname))//'Database'
- write(IMAIN,*) 'make sure file exists'
+ print*,'rank ',myrank,' error opening file: ',prname(1:len_trim(prname))//'Database'
+ print*,'please make sure file exists'
call exit_mpi(myrank,'error opening database file')
endif
!read(IIN,*) nnodes_ext_mesh
@@ -930,7 +930,7 @@
write(IMAIN,*) 'create regions: '
endif
call create_regions_mesh()
-
+
! now done inside create_regions_mesh_ext routine...
! now done inside create_regions_mesh_ext routine...
! Moho boundary parameters, 2-D jacobians and normals
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_coupling_surfaces.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_coupling_surfaces.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -122,7 +122,7 @@
integer :: count_elastic,count_acoustic
! mpi interface communication
- integer, dimension(:), allocatable :: elastic_flag,acoustic_flag !,test_flag
+ integer, dimension(:), allocatable :: elastic_flag,acoustic_flag,test_flag
integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
integer :: max_nibool_interfaces_ext_mesh
logical, dimension(:), allocatable :: mask_ibool
@@ -310,30 +310,34 @@
normal_face(:,i,j) )
! makes sure that it always points away from acoustic element,
! otherwise switch direction
- if( ispec_is_elastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
+ ! note: this should not happen, since we only loop over acoustic elements
+ !if( ispec_is_elastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
+ if( ispec_is_elastic(ispec) ) stop 'error acoustic-elastic coupling surface'
+
enddo
+ enddo
- ! stores informations about this face
- inum = inum + 1
- tmp_ispec(inum) = ispec
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- ! adds all gll points on this face
- igll = igll + 1
+ ! stores informations about this face
+ inum = inum + 1
+ tmp_ispec(inum) = ispec
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds all gll points on this face
+ igll = igll + 1
- ! do we need to store local i,j,k,ispec info? or only global indices iglob?
- tmp_ijk(:,igll,inum) = ijk_face(:,i,j)
+ ! do we need to store local i,j,k,ispec info? or only global indices iglob?
+ tmp_ijk(:,igll,inum) = ijk_face(:,i,j)
- ! stores weighted jacobian and normals
- tmp_jacobian2Dw(igll,inum) = jacobian2Dw_face(i,j)
- tmp_normal(:,igll,inum) = normal_face(:,i,j)
+ ! stores weighted jacobian and normals
+ tmp_jacobian2Dw(igll,inum) = jacobian2Dw_face(i,j)
+ tmp_normal(:,igll,inum) = normal_face(:,i,j)
- ! masks global points ( to avoid redundant counting of faces)
- iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
- mask_ibool(iglob) = .true.
- enddo
+ ! masks global points ( to avoid redundant counting of faces)
+ iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
+ mask_ibool(iglob) = .true.
enddo
+ enddo
! test_flags shouldn't matter, there is only 1 acoustic element touching a coupled surface
! which will be considered in the MPI partition which contains it
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_model.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_model.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -25,8 +25,6 @@
!=====================================================================
-! wrapper function
-
subroutine get_model(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
materials_ext_mesh,nmat_ext_mesh, &
undef_mat_prop,nundefMat_ext_mesh, &
@@ -36,58 +34,7 @@
use create_regions_mesh_ext_par
implicit none
- ! number of spectral elements in each block
- integer :: myrank,nspec
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
- ! external mesh
- integer :: nelmnts_ext_mesh
- integer :: nmat_ext_mesh,nundefMat_ext_mesh
- integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
- character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
- ! anisotropy
- logical :: ANISOTROPY
- character(len=256) LOCAL_PATH
- !----------------------------------------------------------
- ! USER Parameter
- ! daniel: TODO -- uses Piero's get_model_PREM routine rather than default one
- logical,parameter :: USE_PIERO_MODEL = .true.
-
- !----------------------------------------------------------
-
- if( USE_PIERO_MODEL ) then
- ! Using PREM model instead
- call get_model_PREM(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
- materials_ext_mesh,nmat_ext_mesh, &
- undef_mat_prop,nundefMat_ext_mesh, &
- ANISOTROPY,LOCAL_PATH)
-
- else
- ! DEFAULT routine
- call get_model_d(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
- materials_ext_mesh,nmat_ext_mesh, &
- undef_mat_prop,nundefMat_ext_mesh, &
- ANISOTROPY,LOCAL_PATH)
- endif
-
- end subroutine get_model
-
- !
- !-----------------------------------------------------------------------------------------------
- !
-
-! default get model routine
-
- subroutine get_model_d(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
- materials_ext_mesh,nmat_ext_mesh, &
- undef_mat_prop,nundefMat_ext_mesh, &
- ANISOTROPY,LOCAL_PATH)
-
-
- use create_regions_mesh_ext_par
- implicit none
-
! number of spectral elements in each block
integer :: myrank,nspec
@@ -104,6 +51,8 @@
! anisotropy
logical :: ANISOTROPY
+ character(len=256) LOCAL_PATH
+
! local parameters
real(kind=CUSTOM_REAL) :: vp,vs,rho,qmu_atten
real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
@@ -114,24 +63,16 @@
afactor,bfactor,cfactor
integer :: ispec,i,j,k
-
+
! material domain
integer :: idomain_id
-
+
integer :: imaterial_id,imaterial_def
! gll point location
- double precision :: xmesh,ymesh,zmesh
+ double precision :: xmesh,ymesh,zmesh
integer :: iglob
- character(len=256) LOCAL_PATH
- ! variables for importing models from files in SPECFEM format, e.g., proc000000_vp.bin etc.
- ! can be used for importing updated model in iterative inversions
- logical,parameter :: USE_EXTERNAL_FILES = .false.
-
- ! use acoustic domains for simulation
- logical,parameter :: USE_PURE_ACOUSTIC_MOD = .false.
-
! initializes element domain flags
ispec_is_acoustic(:) = .false.
ispec_is_elastic(:) = .false.
@@ -139,7 +80,7 @@
! prepares tomography model if needed for elements with undefined material definitions
if( nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO ) then
- call model_tomography_broadcast(myrank)
+ call model_tomography_broadcast(myrank)
endif
! prepares external model values if needed
@@ -168,7 +109,7 @@
vp = 0._CUSTOM_REAL
vs = 0._CUSTOM_REAL
rho = 0._CUSTOM_REAL
-
+
rho_s = 0._CUSTOM_REAL
kappa_s = 0._CUSTOM_REAL
rho_f = 0._CUSTOM_REAL
@@ -184,7 +125,7 @@
kyy = 0._CUSTOM_REAL
kyz = 0._CUSTOM_REAL
kzz = 0._CUSTOM_REAL
-
+
qmu_atten = 0._CUSTOM_REAL
c11 = 0._CUSTOM_REAL
@@ -216,7 +157,7 @@
zmesh = zstore_dummy(iglob)
! material index 1: associated material number
- ! 1 = acoustic, 2 = elastic, 3 = poroelastic, -1 = undefined tomographic
+ ! 1 = acoustic, 2 = elastic, 3 = poroelastic, -1 = undefined tomographic
imaterial_id = mat_ext_mesh(1,ispec)
! material index 2: associated material definition
@@ -235,12 +176,12 @@
c22,c23,c24,c25,c26,c33, &
c34,c35,c36,c44,c45,c46,c55,c56,c66, &
ANISOTROPY)
-
+
! stores velocity model
- if(idomain_id == IDOMAIN_ACOUSTIC .or. idomain_id == IDOMAIN_ELASTIC) then
-
+ if(idomain_id == IDOMAIN_ACOUSTIC .or. idomain_id == IDOMAIN_ELASTIC) then
+
! elastic or acoustic material
! density
@@ -266,8 +207,8 @@
tortstore(i,j,k,ispec) = 1.d0
!end pll
- else
-
+ else
+
! poroelastic material
! solid properties
@@ -339,11 +280,11 @@
! stores material domain
select case( idomain_id )
- case( IDOMAIN_ACOUSTIC )
+ case( IDOMAIN_ACOUSTIC )
ispec_is_acoustic(ispec) = .true.
case( IDOMAIN_ELASTIC )
ispec_is_elastic(ispec) = .true.
- case( IDOMAIN_POROELASTIC )
+ case( IDOMAIN_POROELASTIC )
ispec_is_poroelastic(ispec) = .true.
case default
stop 'error material domain index'
@@ -386,7 +327,7 @@
if( IMODEL == IMODEL_GLL ) then
! note:
! import the model from files in SPECFEM format
- ! note that those those files should be saved in LOCAL_PATH
+ ! note that those those files should be saved in LOCAL_PATH
call model_gll(myrank,nspec,LOCAL_PATH)
endif
@@ -418,7 +359,7 @@
integer, intent(in) :: nundefMat_ext_mesh
character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
- integer, intent(in) :: imaterial_id,imaterial_def
+ integer, intent(in) :: imaterial_id,imaterial_def
double precision, intent(in) :: xmesh,ymesh,zmesh
@@ -436,14 +377,16 @@
! local parameters
integer :: iflag_aniso
-
+ integer :: iundef,imaterial_PB
+
! use acoustic domains for simulation
logical,parameter :: USE_PURE_ACOUSTIC_MOD = .false.
! initializes with default values
+ ! no anisotropy
iflag_aniso = 0
idomain_id = IDOMAIN_ELASTIC
-
+
! selects chosen velocity model
select case( IMODEL )
@@ -457,374 +400,59 @@
iflag_aniso,qmu_atten,idomain_id, &
rho_s,kappa_s,rho_f,kappa_f,eta_f,kappa_fr,mu_fr, &
phi,tort,kxx,kxy,kxz,kyy,kyz,kzz)
-
+
case( IMODEL_1D_PREM )
! 1D model profile from PREM
call model_1D_prem_iso(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
-
+
+ case( IMODEL_1D_PREM_PB )
+ ! 1D model profile from PREM modified by Piero
+ imaterial_PB = abs(imaterial_id)
+ call model_1D_PREM_routine_PB(xmesh,ymesh,zmesh,rho,vp,vs,imaterial_PB)
+ ! attenuation: arbitrary value, see maximum in constants.h
+ qmu_atten = ATTENUATION_COMP_MAXIMUM
+ ! sets acoustic/elastic domain as given in materials properties
+ iundef = - imaterial_id ! iundef must be positive
+ read(undef_mat_prop(6,iundef),*) idomain_id
+
case( IMODEL_1D_CASCADIA )
! 1D model profile for Cascadia region
call model_1D_cascadia(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
- end subroutine get_model_d
+ case( IMODEL_1D_SOCAL )
+ ! 1D model profile for Southern California
+ call model_1D_socal(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
-!
-!-------------------------------------------------------------------------------------------------
-!
- subroutine get_model_PREM(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
- materials_ext_mesh,nmat_ext_mesh, &
- undef_mat_prop,nundefMat_ext_mesh, &
- ANISOTROPY,LOCAL_PATH)
+ case( IMODEL_SALTON_TROUGH )
+ ! gets model values from tomography file
+ call model_salton_trough(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
- use create_regions_mesh_ext_par
- implicit none
+ case( IMODEL_TOMO )
+ ! gets model values from tomography file
+ call model_tomography(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
- ! number of spectral elements in each block
- integer :: myrank,nspec
+ case( IMODEL_USER_EXTERNAL )
+ ! user model from external routine
+ ! adds/gets velocity model as specified in model_external_values.f90
+ call model_external_values(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten,iflag_aniso,idomain_id)
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ case default
+ stop 'error: model not implemented yet'
+ end select
- ! external mesh
- integer :: nelmnts_ext_mesh
- integer :: nmat_ext_mesh,nundefMat_ext_mesh
+ ! adds anisotropic default model
+ if( ANISOTROPY ) then
+ call model_aniso(iflag_aniso,rho,vp,vs, &
+ c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33, &
+ c34,c35,c36,c44,c45,c46,c55,c56,c66)
+ endif
- integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
- character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
-
- ! anisotropy
- logical :: ANISOTROPY
-
- ! local parameters
- real(kind=CUSTOM_REAL) :: vp,vs,rho,qmu_atten
- real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
- c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
- integer :: ispec,i,j,k,iundef,ier
- integer :: iflag,flag_below,flag_above
- integer :: iflag_aniso,idomain_id,imaterial_id
-
-!PB
- integer :: imaterial_PB
-
- ! gll point location
- double precision :: xloc,yloc,zloc
- integer :: iglob
- character(len=256) LOCAL_PATH,prname_lp
- real, dimension(:,:,:,:),allocatable :: vp_read,vs_read,rho_read
-
- ! variables for importing models from files in SPECFEM format, e.g., proc000000_vp.bin etc.
- ! can be used for importing updated model in iterative inversions
- logical,parameter :: USE_EXTERNAL_FILES = .false.
-
- ! use acoustic domains for simulation
- logical,parameter :: USE_PURE_ACOUSTIC_MOD = .false.
-
-
-
- ! initializes element domain flags
- ispec_is_acoustic(:) = .false.
- ispec_is_elastic(:) = .false.
- ispec_is_poroelastic(:) = .false.
-
-!PB AT PRESENT I DON'T NEED THIS CALL CAUSE I'M USING PREM MODEL
-
- ! prepares tomography model if needed for elements with undefined material definitions
-! if( nundefMat_ext_mesh > 0 ) then
-! call model_tomography_broadcast(myrank)
-! endif
-
- ! prepares external model values if needed
- if( USE_MODEL_EXTERNAL_VALUES ) then
- call model_external_broadcast(myrank)
+ ! for pure acoustic simulations (a way of avoiding re-mesh, re-partition etc.)
+ ! can be used to compare elastic & acoustic reflections in exploration seismology
+ ! do NOT use it unless you are confident
+ if( USE_PURE_ACOUSTIC_MOD ) then
+ idomain_id = IDOMAIN_ACOUSTIC
endif
-! ! Piero, read bedrock file
-! in case, see file model_interface_bedrock.f90:
-! call model_bedrock_broadcast(myrank)
-
-
- ! material properties on all GLL points: taken from material values defined for
- ! each spectral element in input mesh
- do ispec = 1, nspec
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
-
- ! material index 1: associated material number
- imaterial_id = mat_ext_mesh(1,ispec)
-
- ! check if the material is known or unknown
- if( imaterial_id > 0) then
- ! gets velocity model as specified by (cubit) mesh files
-
- ! density
- ! materials_ext_mesh format:
- ! #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
- rho = materials_ext_mesh(1,imaterial_id)
-
- ! isotropic values: vp, vs
- vp = materials_ext_mesh(2,imaterial_id)
- vs = materials_ext_mesh(3,imaterial_id)
-
- ! attenuation
- qmu_atten = materials_ext_mesh(4,imaterial_id)
-
- ! anisotropy
- iflag_aniso = materials_ext_mesh(5,imaterial_id)
-
- ! material domain_id
- idomain_id = materials_ext_mesh(6,imaterial_id)
-
- else if (mat_ext_mesh(2,ispec) == 1) then
-
- stop 'material: interface not implemented yet'
-
- do iundef = 1,nundefMat_ext_mesh
- if(trim(undef_mat_prop(2,iundef)) == 'interface') then
- read(undef_mat_prop(3,iundef),'(1i3)') flag_below
- read(undef_mat_prop(4,iundef),'(1i3)') flag_above
- endif
- enddo
-
- ! see file model_interface_bedrock.f90: routine interface()
- !call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
-
- ! dummy: takes 1. defined material
- iflag = 1
- rho = materials_ext_mesh(1,iflag)
- vp = materials_ext_mesh(2,iflag)
- vs = materials_ext_mesh(3,iflag)
- qmu_atten = materials_ext_mesh(4,iflag)
- iflag_aniso = materials_ext_mesh(5,iflag)
- idomain_id = materials_ext_mesh(6,iflag)
-
- else if ( mat_ext_mesh(2,ispec) == 2 ) then
-
- imaterial_PB = abs(imaterial_id)
- ! material definition undefined, uses definition from tomography model
- ! GLL point location
- iglob = ibool(i,j,k,ispec)
- xloc = xstore_dummy(iglob)
- yloc = ystore_dummy(iglob)
- zloc = zstore_dummy(iglob)
-
- call PREM_routine(xloc,yloc,zloc, &
- rho,vp,vs,imaterial_PB)
-
-
-!PB COMMENTED THE CALL TO
-!model_tomography SINCE I WANT
-!EACH GLL POINT HAVING
-!PB VELOCITY VALUE DEFINED BY THE
-!PREM_routine
- ! gets model values from tomography file
-! call model_tomography(xloc,yloc,zloc, &
-! rho,vp,vs)
-
- qmu_atten = ATTENUATION_COMP_MAXIMUM ! attenuation: arbitrary value, see maximum in constants.h
- iflag_aniso = 0 ! no anisotropy
-
- ! sets acoustic/elastic domain as given in materials properties
- iundef = - imaterial_id ! iundef must be positive
- read(undef_mat_prop(6,iundef),*) idomain_id
- ! or
- !idomain_id = IDOMAIN_ELASTIC ! forces to be elastic domain
-
-!PB
-
- else
-
- stop 'material: not implemented yet'
-
- end if
-
- ! adds/gets velocity model as specified in model_external_values.f90
- if( USE_MODEL_EXTERNAL_VALUES ) then
- call model_external_values(i,j,k,ispec,idomain_id,imaterial_id, &
- nspec,ibool, &
- iflag_aniso,qmu_atten, &
- rho,vp,vs, &
- c11,c12,c13,c14,c15,c16, &
- c22,c23,c24,c25,c26,c33, &
- c34,c35,c36,c44,c45,c46, &
- c55,c56,c66,ANISOTROPY)
- endif
-
- ! adds anisotropic default model
- if( ANISOTROPY .and. .not. USE_MODEL_EXTERNAL_VALUES ) then
- call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
- c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45, &
- c46,c55,c56,c66)
-
- endif
-
- ! stores velocity model
-
- ! density
- rhostore(i,j,k,ispec) = rho
-
- ! kappa, mu
- kappastore(i,j,k,ispec) = rho*( vp*vp - FOUR_THIRDS*vs*vs )
- mustore(i,j,k,ispec) = rho*vs*vs
-
- ! attenuation
- qmu_attenuation_store(i,j,k,ispec) = qmu_atten
-
- ! Stacey, a completer par la suite
- rho_vp(i,j,k,ispec) = rho*vp
- rho_vs(i,j,k,ispec) = rho*vs
- !end pll
-
- ! adds anisotropic perturbation to vp, vs
- if( ANISOTROPY ) then
- !call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
- ! c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
- c11store(i,j,k,ispec) = c11
- c12store(i,j,k,ispec) = c12
- c13store(i,j,k,ispec) = c13
- c14store(i,j,k,ispec) = c14
- c15store(i,j,k,ispec) = c15
- c16store(i,j,k,ispec) = c16
- c22store(i,j,k,ispec) = c22
- c23store(i,j,k,ispec) = c23
- c24store(i,j,k,ispec) = c24
- c25store(i,j,k,ispec) = c25
- c26store(i,j,k,ispec) = c26
- c33store(i,j,k,ispec) = c33
- c34store(i,j,k,ispec) = c34
- c35store(i,j,k,ispec) = c35
- c36store(i,j,k,ispec) = c36
- c44store(i,j,k,ispec) = c44
- c45store(i,j,k,ispec) = c45
- c46store(i,j,k,ispec) = c46
- c55store(i,j,k,ispec) = c55
- c56store(i,j,k,ispec) = c56
- c66store(i,j,k,ispec) = c66
- endif
-
- ! for pure acoustic simulations (a way of avoiding re-mesh, re-partition etc.)
- ! can be used to compare elastic & acoustic reflections in exploration seismology
- ! do NOT use it unless you are confident
- if( USE_PURE_ACOUSTIC_MOD ) then
- idomain_id = IDOMAIN_ACOUSTIC
- endif
-
- ! material domain
- !print*,'velocity model:',ispec,idomain_id
- if( idomain_id == IDOMAIN_ACOUSTIC ) then
- ispec_is_acoustic(ispec) = .true.
- else if( idomain_id == IDOMAIN_ELASTIC ) then
- ispec_is_elastic(ispec) = .true.
- else if( idomain_id == IDOMAIN_POROELASTIC ) then
- stop 'poroelastic material domain not implemented yet'
- ispec_is_poroelastic(ispec) = .true.
- else
- stop 'error material domain index'
- endif
-
- enddo
- enddo
- enddo
- !print*,myrank,'ispec:',ispec,'rho:',rhostore(1,1,1,ispec),'vp:',vpstore(1,1,1,ispec),'vs:',vsstore(1,1,1,ispec)
- enddo
-
- ! checks material domains
- do ispec=1,nspec
- ! checks if domain is set
- if( (ispec_is_acoustic(ispec) .eqv. .false.) &
- .and. (ispec_is_elastic(ispec) .eqv. .false.) &
- .and. (ispec_is_poroelastic(ispec) .eqv. .false.) ) then
- print*,'error material domain not assigned to element:',ispec
- print*,'acoustic: ',ispec_is_acoustic(ispec)
- print*,'elastic: ',ispec_is_elastic(ispec)
- print*,'poroelastic: ',ispec_is_poroelastic(ispec)
- stop 'error material domain index element'
- endif
- ! checks if domain is unique
- if( ispec_is_acoustic(ispec) .eqv. .true. .and. ispec_is_elastic(ispec) .eqv. .true. ) then
- print*,'error material domain assigned twice to element:',ispec
- print*,'acoustic: ',ispec_is_acoustic(ispec)
- print*,'elastic: ',ispec_is_elastic(ispec)
- print*,'poroelastic: ',ispec_is_poroelastic(ispec)
- stop 'error material domain index element'
- endif
- enddo
-
-! !! DK DK store the position of the six stations to be able to
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! in case, see file model_interface_bedrock.f90: routine model_bedrock_store()
-
-
-! import the model from files in SPECFEM format
-! note that those those files should be saved in LOCAL_PATH
-
- if( USE_EXTERNAL_FILES ) then
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!! if only vp structure is available (as is often the case in exploration seismology),
-!!! use lines for vp only
-
- ! processors name
- write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
-
- allocate( rho_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
- if( ier /= 0 ) stop 'error allocating array rho_read'
- write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
- open(unit=28,file=prname_lp(1:len_trim(prname_lp))//'rho.bin',&
- status='unknown',action='read',form='unformatted')
- read(28) rho_read
- close(28)
-
- allocate( vp_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
- if( ier /= 0 ) stop 'error allocating array vp_read'
- write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
- open(unit=28,file=prname_lp(1:len_trim(prname_lp))//'vp.bin',&
- status='unknown',action='read',form='unformatted')
- read(28) vp_read
- close(28)
-
- allocate( vs_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
- if( ier /= 0 ) stop 'error allocating array vs_read'
- write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
- open(unit=28,file=prname_lp(1:len_trim(prname_lp))//'vs.bin',&
- status='unknown',action='read',form='unformatted')
- read(28) vs_read
- close(28)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!! in cases where density structure is not given
-!!! modify according to your desire
-
-! rho_read = 1000.0
-! where ( mustore > 100.0 ) &
-! rho_read = (1.6612 * (vp_read / 1000.0) &
-! -0.4720 * (vp_read / 1000.0)**2 &
-! +0.0671 * (vp_read / 1000.0)**3 &
-! -0.0043 * (vp_read / 1000.0)**4 &
-! +0.000106*(vp_read / 1000.0)**5 )*1000.0
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!! in cases where shear wavespeed structure is not given
-!!! modify according to your desire
-
-! vs_read = 0.0
-! where ( mustore > 100.0 ) vs_read = vp_read / sqrt(3.0)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!! update arrays that will be saved and used in the solver xspecfem3D
-!!! the following part is neccessary if you uncommented something above
-
- rhostore = rho_read
- kappastore = rhostore * ( vp_read * vp_read - FOUR_THIRDS * vs_read * vs_read )
- mustore = rhostore * vs_read * vs_read
- rho_vp = rhostore * vp_read
- rho_vs = rhostore * vs_read
-
- ! free memory
- deallocate( rho_read,vp_read,vs_read)
-
- endif ! USE_EXTERNAL_FILES
-
- end subroutine get_model_PREM
-
+ end subroutine get_model_values
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -683,6 +683,7 @@
! loop on all the colors to determine the color with the smallest number
! of elements and for which there is no conflict
nb_elems_in_color_chosen = 2147000000 ! start with extremely large unrealistic value
+ icolor_chosen = 0
do icolor = icolormin,icolormax
if (.not. icolor_conflict_found(icolor) .and. nb_elems_in_this_color(icolor) < nb_elems_in_color_chosen) then
icolor_chosen = icolor
@@ -868,6 +869,7 @@
! loop on all the colors to determine the color with the smallest number of elements
! and for which there is no conflict
nb_elems_in_color_chosen = 2147000000 ! start with extremely large unrealistic value
+ icolor_chosen = 0
do icolor_target = icolormin,icolormax
if (.not. icolor_conflict_found(icolor_target) .and. &
nb_elems_in_this_color(icolor_target) < nb_elems_in_color_chosen) then
@@ -995,6 +997,7 @@
! loops on all the colors to determine the color with the smallest number of elements
! and for which there is no conflict
nb_elems_in_color_chosen = 2147000000 ! start with extremely large unrealistic value
+ icolor_chosen = 0
do icolor_target = icolormin,icolormax
if (.not. icolor_conflict_found(icolor_target) .and. &
nb_elems_in_this_color(icolor_target) < nb_elems_in_color_chosen) then
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/memory_eval.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/memory_eval.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -166,7 +166,7 @@
max_interface_size_ext_mesh,nspec2D_xmin,nspec2D_xmax, &
nspec2D_ymin,nspec2D_ymax,nspec2D_bottom,nspec2D_top
- integer,intent(inout) :: static_memory_size_request
+ double precision,intent(inout) :: static_memory_size_request
! local parameters
integer :: static_memory_size
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_cascadia.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_cascadia.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_cascadia.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -48,7 +48,7 @@
! local parameters
real(kind=CUSTOM_REAL) :: x,y,z
real(kind=CUSTOM_REAL) :: depth
- real(kind=CUSTOM_REAL) :: elevation,distmin
+ real(kind=CUSTOM_REAL) :: elevation,distmin
! converts GLL point location to real
x = xmesh
@@ -61,19 +61,19 @@
call get_topo_elevation_free_closest(x,y,elevation,distmin, &
nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
-
+
! depth in Z-direction
- if( distmin < HUGEVAL ) then
+ if( distmin < HUGEVAL ) then
depth = elevation - z
else
depth = - z
endif
! depth in km
- depth = depth / 1000.0
+ depth = depth / 1000.0
! 1D profile Cascadia
-
+
! super-imposes values
if( depth < 1.0 ) then
! vp in m/s
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_prem.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_prem.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_prem.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -51,29 +51,29 @@
implicit none
! GLL point location
- double precision, intent(in) :: xmesh,ymesh,zmesh
-
+ double precision, intent(in) :: xmesh,ymesh,zmesh
+
! material properties
real(kind=CUSTOM_REAL), intent(inout) :: rho_prem,vp_prem,vs_prem,qmu_atten
-
+
! local parameters
real(kind=CUSTOM_REAL) :: xloc,yloc,zloc
real(kind=CUSTOM_REAL) :: depth
- real(kind=CUSTOM_REAL) :: elevation,distmin
+ real(kind=CUSTOM_REAL) :: elevation,distmin
double precision :: x,rho,drhodr,vp,vs,Qkappa,Qmu
double precision :: &
- R_EARTH,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R_EARTH_M,RICB,RCMB,RTOPDDOUBLEPRIME, &
R771,R600,R670,R400,R220,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
double precision :: r
! uses crustal values from other models (like crust2.0) than prem
! set to .false. to use PREM crustal values, otherwise will take mantle values up to surface
- logical,parameter :: CRUSTAL = .false.
+ logical,parameter :: CRUSTAL = .false.
! avoids crustal values, uses Moho values up to the surface
logical,parameter :: SUPPRESS_CRUSTAL_MESH = .false.
! same properties everywhere in PREM crust if we decide to define only one layer in the crust
logical,parameter :: ONE_CRUST = .false.
-
+
! GLL point location converted to real
xloc = xmesh
yloc = ymesh
@@ -85,16 +85,16 @@
call get_topo_elevation_free_closest(xloc,yloc,elevation,distmin, &
nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
-
+
! depth in Z-direction
- if( distmin < HUGEVAL ) then
+ if( distmin < HUGEVAL ) then
depth = elevation - zloc
else
depth = - zloc
endif
! PREM layers (in m)
- R_EARTH = 6371000.d0
+ R_EARTH_M = 6371000.d0
ROCEAN = 6368000.d0
RMIDDLE_CRUST = 6356000.d0
RMOHO = 6346600.d0
@@ -113,10 +113,18 @@
! normalized radius
x = r / R_EARTH
-
+
! given a normalized radius x, gives the non-dimensionalized density rho,
! speeds vp and vs, and the quality factors Qkappa and Qmu
+ ! initializes
+ rho = 0.d0
+ vp = 0.d0
+ vs = 0.d0
+ Qmu = 0.d0
+ Qkappa = 0.d0
+ drhodr = 0.d0
+
!
!--- inner core
!
@@ -271,5 +279,102 @@
vs_prem=vs*1000.0d0
qmu_atten = Qmu
-
+
end subroutine model_1D_prem_iso
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!PB ROUTINE WHICH ASSINGS AT EACH GLL POINT THE VALUES OF v_p v_s AND rho TAKEN FROM THE PREM MODEL
+!PB IT'S CALLED IN get_model.f90(:149) WITHIN THE generate_databases PROCESS
+
+ subroutine model_1D_PREM_routine_PB(xloc,yloc,zloc,ro_prem,vp_prem,vs_prem,idom)
+
+ use generate_databases_par,only: CUSTOM_REAL
+ implicit none
+
+ double precision, intent(in) :: xloc,yloc,zloc
+ real(kind=CUSTOM_REAL),intent(inout) :: ro_prem,vp_prem,vs_prem
+ integer, intent(in) :: idom
+
+ ! local parameters
+ double precision :: r0,r,x_prem
+ !double precision :: ro_prem,vp_prem,vs_prem
+ !character(len=3), intent(in) :: param !rho, vs,vp
+
+ r0=sqrt(xloc**2+yloc**2+zloc**2)
+ ! print*,'xloc,yloc,zloc,r0,idom',xloc,yloc,zloc,r0,idom
+ r=r0/1000.
+
+ x_prem=r/6371. ! Radius (normalized to x(surface)=1 )
+ ! print*,'x_prem',x_prem
+ IF(idom==1)THEN ! upper crustal layer
+ !print*,'I am in domain 1'
+ ro_prem=2.6
+ vp_prem=5.8
+ vs_prem=3.2
+ ! print*,'ro,vp,vs form domain 1',ro_prem,vp_prem,vs_prem
+ ELSEIF(idom==2)THEN
+ ro_prem=2.9 ! lower crustal layer
+ vp_prem=6.8
+ vs_prem=3.9
+ ELSEIF(idom==3)THEN
+ ro_prem=2.691+.6924*x_prem ! upper mantle
+ vp_prem=4.1875+3.9382*x_prem
+ vs_prem=2.1519+2.3481*x_prem
+ ELSEIF(idom==4)THEN
+ ro_prem=7.1089-3.8045*x_prem
+ vp_prem=20.3926-12.2569*x_prem
+ vs_prem=8.9496-4.4597*x_prem
+ ELSEIF(idom==5)THEN
+ ro_prem=11.2494-8.0298*x_prem
+ vp_prem=39.7027-32.6166*x_prem
+ vs_prem=22.3512-18.5856*x_prem
+ ELSEIF(idom==6)THEN
+ ro_prem=5.3197-1.4836*x_prem
+ vp_prem=19.0957-9.8672*x_prem
+ vs_prem=9.9839-4.9324*x_prem
+ ELSEIF(idom==7)THEN !lower mantle
+ ro_prem=7.9565-6.4761*x_prem+5.5283*x_prem**2-3.0807*x_prem**3
+ vp_prem=29.2766-23.6027*x_prem+5.5242*x_prem**2-2.5514*x_prem**3
+ vs_prem=22.3459-17.2473*x_prem-2.0834*x_prem**2+0.9783*x_prem**3
+ ELSEIF(idom==8)THEN
+ ro_prem=7.9565-6.4761*x_prem+5.5283*x_prem**2-3.0807*x_prem**3
+ vp_prem=24.9520-40.4673*x_prem+51.4832*x_prem**2-26.6419*x_prem**3
+ vs_prem=11.1671-13.7818*x_prem+17.4575*x_prem**2-9.2777*x_prem**3
+ ELSEIF(idom==9)THEN
+ ro_prem=7.9565-6.4761*x_prem+5.5283*x_prem**2-3.0807*x_prem**3
+ vp_prem=15.3891-5.3181*x_prem+5.5242*x_prem**2-2.5514*x_prem**3
+ vs_prem=6.9254+1.4672*x_prem-2.0834*x_prem**2+.9783*x_prem**3
+ ELSEIF(idom==10)THEN ! outer core
+ ro_prem=12.5815-1.2638*x_prem-3.6426*x_prem**2-5.5281*x_prem**3
+ vp_prem=11.0487-4.0362*x_prem+4.8023*x_prem**2-13.5732*x_prem**3
+ vs_prem=0.00
+ ELSEIF(idom==11)THEN ! inner core
+ ro_prem=13.0885-8.8381*x_prem**2
+ vp_prem=11.2622-6.3640*x_prem**2
+ vs_prem=3.6678-4.4475*x_prem**2
+ ENDIF
+
+ ! print*,'ro,vp,vs passed from routine and not multiplied',ro_prem,vp_prem,vs_prem
+
+ ro_prem=ro_prem*1000
+ vp_prem=vp_prem*1000
+ vs_prem=vs_prem*1000
+
+ ! print*,'ro,vp,vs passed from PREM_ROUTINE multiplied by 1000',ro_prem,vp_prem,vs_prem
+
+ ! if (param=='rho') then
+ ! prem_sub=ro_prem*1000.
+ ! elseif (param=='v_p') then
+ ! prem_sub=vp_prem*1000.
+ ! elseif (param=='v_s') then
+ ! prem_sub=vs_prem*1000.
+ ! else
+ ! write(6,*)'ERROR IN PREM_SUB FUNCTION:',param,'NOT AN OPTION'
+ ! stop
+
+ end subroutine model_1D_PREM_routine_PB
+
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_socal.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_socal.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_1d_socal.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -26,13 +26,13 @@
!--------------------------------------------------------------------------------------------------
!
-! 1D Southern California model
+! 1D Southern California model
!
! model is the standard model used in southern California:
-! Kanamori and Hadley (1975), Dreger and Helmberger (1990), Wald-Hutton,Given (1995)
+! Kanamori and Hadley (1975), Dreger and Helmberger (1990), Wald-Hutton,Given (1995)
!
!--------------------------------------------------------------------------------------------------
-
+
subroutine model_1D_socal(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
! given a GLL point, returns super-imposed velocity model values
@@ -49,12 +49,12 @@
! local parameters
real(kind=CUSTOM_REAL) :: depth,x,y,z
-
+
! mesh point location
x = xmesh
y = ymesh
z = zmesh
-
+
! depth in m
depth = -zmesh
@@ -88,5 +88,5 @@
! no attenuation information
qmu_atten = 0.d0
-
+
end subroutine model_1D_socal
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_default.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_default.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_default.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -49,7 +49,7 @@
integer, intent(in) :: nundefMat_ext_mesh
character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
- integer, intent(in) :: imaterial_id,imaterial_def
+ integer, intent(in) :: imaterial_id,imaterial_def
double precision, intent(in) :: xmesh,ymesh,zmesh
@@ -60,11 +60,11 @@
real(kind=CUSTOM_REAL) :: kappa_s,kappa_f,kappa_fr,mu_fr,rho_s,rho_f,phi,tort,eta_f, &
kxx,kxy,kxz,kyy,kyz,kzz
-
+
! local parameters
integer :: iflag,flag_below,flag_above
integer :: iundef
-
+
! check if the material is known or unknown
if( imaterial_id > 0 ) then
! gets velocity model as specified by (cubit) mesh files for elastic & acoustic
@@ -74,8 +74,8 @@
idomain_id = materials_ext_mesh(6,imaterial_id)
select case( idomain_id )
-
- case( IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC)
+
+ case( IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC)
! elastic or acoustic
! density
@@ -93,7 +93,7 @@
! anisotropy
iflag_aniso = materials_ext_mesh(5,imaterial_id)
- case( IDOMAIN_POROELASTIC )
+ case( IDOMAIN_POROELASTIC )
! poroelastic
! materials_ext_mesh format:
! rho_s,kappa_s,rho_f,kappa_f,eta_f,kappa_fr,mu_fr,phi,tort,kxx,kxy,kxz,kyy,kyz,kzz
@@ -120,9 +120,9 @@
case default
print*,'error: domain id = ',idomain_id,'not recognized'
stop 'error: domain not recognized'
-
+
end select
-
+
else if ( imaterial_def == 1 ) then
stop 'material: interface not implemented yet'
@@ -142,7 +142,7 @@
rho = materials_ext_mesh(1,iflag)
vp = materials_ext_mesh(2,iflag)
vs = materials_ext_mesh(3,iflag)
- qmu_atten = materials_ext_mesh(4,iflag)
+ qmu_atten = materials_ext_mesh(4,iflag)
iflag_aniso = materials_ext_mesh(5,iflag)
idomain_id = materials_ext_mesh(6,iflag)
@@ -154,7 +154,7 @@
call model_tomography(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
! no anisotropy
- iflag_aniso = 0
+ iflag_aniso = 0
! sets acoustic/elastic domain as given in materials properties
iundef = - imaterial_id ! iundef must be positive
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_external_values.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_external_values.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_external_values.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -124,7 +124,7 @@
use create_regions_mesh_ext_par
implicit none
- ! GLL point
+ ! GLL point
double precision, intent(in) :: xmesh,ymesh,zmesh
! density, Vp and Vs
@@ -168,14 +168,14 @@
call get_topo_elevation_free_closest(x,y,elevation,distmin, &
nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
-
+
! depth in Z-direction
- if( distmin < HUGEVAL ) then
+ if( distmin < HUGEVAL ) then
depth = elevation - z
else
depth = zmax - z
endif
-
+
! normalizes depth between 0 and 1
if( abs( zmax - zmin ) > TINYVAL ) depth = depth / (zmax - zmin)
@@ -200,6 +200,6 @@
iflag_aniso = 0
! elastic material
- idomain_id = IDOMAIN_ELASTIC
-
+ idomain_id = IDOMAIN_ELASTIC
+
end subroutine model_external_values
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_gll.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_gll.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_gll.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -39,14 +39,14 @@
use create_regions_mesh_ext_par
implicit none
-
+
integer, intent(in) :: myrank,nspec
character(len=256) :: LOCAL_PATH
-
- ! local parameters
+
+ ! local parameters
real, dimension(:,:,:,:),allocatable :: vp_read,vs_read,rho_read
integer :: ier
- character(len=256) :: prname_lp,filename
+ character(len=256) :: prname_lp,filename
! processors name
write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
@@ -81,7 +81,7 @@
filename = prname_lp(1:len_trim(prname_lp))//'vs.bin'
open(unit=28,file=trim(filename),status='unknown',action='read',form='unformatted')
-
+
read(28) vs_read
close(28)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_salton_trough.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_salton_trough.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_salton_trough.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -43,7 +43,7 @@
double precision, parameter :: GOCAD_ST_V_X = 109670.74, GOCAD_ST_V_Y = 71530.72
double precision, parameter :: GOCAD_ST_W_Z = 7666.334
double precision, parameter :: GOCAD_ST_NO_DATA_VALUE = -99999
-
+
real,dimension(:,:,:),allocatable :: vp_array
end module salton_trough_par
@@ -63,8 +63,8 @@
! local parameters
integer :: ier
-
-
+
+
allocate(vp_array(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating vp_array for salton')
@@ -92,11 +92,11 @@
! array length
reclen=(GOCAD_ST_NU * GOCAD_ST_NV * GOCAD_ST_NW) * 4
-
- ! file name
+
+ ! file name
call get_value_string(SALTON_SEA_MODEL_FILE,'model.SALTON_SEA_MODEL_FILE', &
'DATA/st_3D_block_harvard/regrid3_vel_p.bin')
-
+
! reads in file values
open(11,file=trim(SALTON_SEA_MODEL_FILE), &
status='old',action='read',form='unformatted',access='direct',recl=reclen,iostat=ios)
@@ -104,10 +104,10 @@
print *,'error opening file: ',trim(SALTON_SEA_MODEL_FILE),' iostat = ', ios
call exit_mpi(0,'Error opening file salton trough')
endif
-
- read(11,rec=1,iostat=ios) vp_array
+
+ read(11,rec=1,iostat=ios) vp_array
if (ios /= 0) stop 'Error reading vp_array'
-
+
close(11)
end subroutine read_salton_sea_model
@@ -126,7 +126,7 @@
use create_regions_mesh_ext_par
implicit none
- ! GLL point
+ ! GLL point
double precision, intent(in) :: xmesh,ymesh,zmesh
! density, Vp and Vs
@@ -139,20 +139,20 @@
double precision :: uc,vc,wc
double precision :: vp_st,vs_st,rho_st
- ! GLL point location converted to u,v,w
+ ! GLL point location converted to u,v,w
call vx_xyz2uvw(xmesh,ymesh,zmesh,uc,vc,wc)
! model values
call vx_xyz_interp(uc,vc,wc,vp_st,vs_st,rho_st)
-
+
! converts to custom real
vp = vp_st
vs = vs_st
rho = rho_st
-
+
! no attenuation info
qmu_atten = 0.0
-
+
end subroutine model_salton_trough
!
@@ -221,7 +221,7 @@
v8 = vp_array(i,j+1,k+1)
vi = vp_array(i+ixi,j+ieta,k+iga)
! print *, v1, v2, v3, v4, v5, v6, v7, v8
-
+
if ((v1 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
(v2 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
(v3 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
@@ -240,7 +240,7 @@
v8 * (1-xi) * eta * ga)
else if ((vi - GOCAD_ST_NO_DATA_VALUE) > eps) then
vp = dble(vi)
-
+
! else if ((v1 - GOCAD_ST_NO_DATA_VALUE) > eps) then
! vp = dble(v1)
! else if ((v2 - GOCAD_ST_NO_DATA_VALUE) > eps) then
@@ -257,21 +257,21 @@
! vp = dble(v7)
! else if ((v7 - GOCAD_ST_NO_DATA_VALUE) > eps) then
! vp = dble(v8)
-
+
else
vp = GOCAD_ST_NO_DATA_VALUE
endif
-
+
! depth
zmesh = wc / (GOCAD_ST_NW - 1) * GOCAD_ST_W_Z + GOCAD_ST_O_Z
-
+
! vs
if (zmesh > -8500.) then
vs = vp / (2 - (0.27*zmesh/(-8500)))
else
vs = vp/1.73
endif
-
+
! density
if (vp > 2160.) then
rho = vp/3 + 1280.
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_tomography.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/model_tomography.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -82,7 +82,7 @@
!if(myrank == 0) call read_external_model()
! broadcast the information read on the master to the nodes, e.g.
!call bcast_all_i(nrecord,1)
-
+
!if( myrank /= 0 ) then
! allocate( vp_tomography(1:nrecord) ,stat=ier)
! if( ier /= 0 ) stop 'error allocating array vp_tomography'
@@ -161,101 +161,9 @@
end subroutine read_model_tomography
-
!
!-------------------------------------------------------------------------------------------------
-!PB ROUTINE WHICH ASSINGS AT EACH GLL POINT THE VALUES OF v_p v_s AND rho TAKEN FROM THE PREM MODEL
-!PB IT'S CALLED IN get_model.f90(:149) WITHIN THE generate_databases PROCESS
-
- subroutine PREM_routine(xloc,yloc,zloc,ro_prem,vp_prem,vs_prem,idom)
-
-double precision, intent(in) :: xloc,yloc,zloc
-integer, intent(in) :: idom
-double precision :: r0,r,x_prem
-!double precision :: ro_prem,vp_prem,vs_prem
-real :: ro_prem,vp_prem,vs_prem
-!character(len=3), intent(in) :: param !rho, vs,vp
-
- r0=sqrt(xloc**2+yloc**2+zloc**2)
-! print*,'xloc,yloc,zloc,r0,idom',xloc,yloc,zloc,r0,idom
- r=r0/1000.
-
- x_prem=r/6371. ! Radius (normalized to x(surface)=1 )
-! print*,'x_prem',x_prem
- IF(idom==1)THEN ! upper crustal layer
- !print*,'I am in domain 1'
- ro_prem=2.6
- vp_prem=5.8
- vs_prem=3.2
-! print*,'ro,vp,vs form domain 1',ro_prem,vp_prem,vs_prem
- ELSEIF(idom==2)THEN
- ro_prem=2.9 ! lower crustal layer
- vp_prem=6.8
- vs_prem=3.9
- ELSEIF(idom==3)THEN
- ro_prem=2.691+.6924*x_prem ! upper mantle
- vp_prem=4.1875+3.9382*x_prem
- vs_prem=2.1519+2.3481*x_prem
- ELSEIF(idom==4)THEN
- ro_prem=7.1089-3.8045*x_prem
- vp_prem=20.3926-12.2569*x_prem
- vs_prem=8.9496-4.4597*x_prem
- ELSEIF(idom==5)THEN
- ro_prem=11.2494-8.0298*x_prem
- vp_prem=39.7027-32.6166*x_prem
- vs_prem=22.3512-18.5856*x_prem
- ELSEIF(idom==6)THEN
- ro_prem=5.3197-1.4836*x_prem
- vp_prem=19.0957-9.8672*x_prem
- vs_prem=9.9839-4.9324*x_prem
- ELSEIF(idom==7)THEN !lower mantle
- ro_prem=7.9565-6.4761*x_prem+5.5283*x_prem**2-3.0807*x_prem**3
- vp_prem=29.2766-23.6027*x_prem+5.5242*x_prem**2-2.5514*x_prem**3
- vs_prem=22.3459-17.2473*x_prem-2.0834*x_prem**2+0.9783*x_prem**3
- ELSEIF(idom==8)THEN
- ro_prem=7.9565-6.4761*x_prem+5.5283*x_prem**2-3.0807*x_prem**3
- vp_prem=24.9520-40.4673*x_prem+51.4832*x_prem**2-26.6419*x_prem**3
- vs_prem=11.1671-13.7818*x_prem+17.4575*x_prem**2-9.2777*x_prem**3
- ELSEIF(idom==9)THEN
- ro_prem=7.9565-6.4761*x_prem+5.5283*x_prem**2-3.0807*x_prem**3
- vp_prem=15.3891-5.3181*x_prem+5.5242*x_prem**2-2.5514*x_prem**3
- vs_prem=6.9254+1.4672*x_prem-2.0834*x_prem**2+.9783*x_prem**3
- ELSEIF(idom==10)THEN ! outer core
- ro_prem=12.5815-1.2638*x_prem-3.6426*x_prem**2-5.5281*x_prem**3
- vp_prem=11.0487-4.0362*x_prem+4.8023*x_prem**2-13.5732*x_prem**3
- vs_prem=0.00
- ELSEIF(idom==11)THEN ! inner core
- ro_prem=13.0885-8.8381*x_prem**2
- vp_prem=11.2622-6.3640*x_prem**2
- vs_prem=3.6678-4.4475*x_prem**2
- ENDIF
-
-! print*,'ro,vp,vs passed from routine and not multiplied',ro_prem,vp_prem,vs_prem
-
- ro_prem=ro_prem*1000
- vp_prem=vp_prem*1000
- vs_prem=vs_prem*1000
-
-! print*,'ro,vp,vs passed from PREM_ROUTINE multiplied by 1000',ro_prem,vp_prem,vs_prem
-
-! if (param=='rho') then
-! prem_sub=ro_prem*1000.
-! elseif (param=='v_p') then
-! prem_sub=vp_prem*1000.
-! elseif (param=='v_s') then
-! prem_sub=vs_prem*1000.
-! else
-! write(6,*)'ERROR IN PREM_SUB FUNCTION:',param,'NOT AN OPTION'
-! stop
-
-end subroutine PREM_routine
-
-
-
-
!
-!-------------------------------------------------------------------------------------------------
-!
subroutine model_tomography(xmesh,ymesh,zmesh,rho_final,vp_final,vs_final,qmu_atten)
@@ -270,7 +178,7 @@
!double precision, intent(in) :: VP_MIN,VS_MIN,RHO_MIN,VP_MAX,VS_MAX,RHO_MAX
double precision, intent(in) :: xmesh,ymesh,zmesh
-
+
real(kind=CUSTOM_REAL), intent(out) :: vp_final,vs_final,rho_final,qmu_atten
! local parameters
@@ -458,6 +366,6 @@
if(rho_final > RHO_MAX) rho_final = RHO_MAX
! attenuation: arbitrary value, see maximum in constants.h
- qmu_atten = ATTENUATION_COMP_MAXIMUM
+ qmu_atten = ATTENUATION_COMP_MAXIMUM
end subroutine model_tomography
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/save_arrays_solver.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/save_arrays_solver.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -66,7 +66,11 @@
nspec_outer_acoustic,nspec_outer_elastic,nspec_outer_poroelastic, &
num_phase_ispec_acoustic,phase_ispec_inner_acoustic, &
num_phase_ispec_elastic,phase_ispec_inner_elastic, &
- num_phase_ispec_poroelastic,phase_ispec_inner_poroelastic)
+ num_phase_ispec_poroelastic,phase_ispec_inner_poroelastic, &
+ num_colors_outer_acoustic,num_colors_inner_acoustic, &
+ num_elem_colors_acoustic, &
+ num_colors_outer_elastic,num_colors_inner_elastic, &
+ num_elem_colors_elastic)
implicit none
@@ -176,6 +180,14 @@
integer :: num_phase_ispec_poroelastic
integer,dimension(num_phase_ispec_poroelastic,2) :: phase_ispec_inner_poroelastic
+ ! mesh coloring
+ integer :: num_colors_outer_acoustic,num_colors_inner_acoustic
+ integer, dimension(num_colors_outer_acoustic + num_colors_inner_acoustic) :: &
+ num_elem_colors_acoustic
+ integer :: num_colors_outer_elastic,num_colors_inner_elastic
+ integer, dimension(num_colors_outer_elastic + num_colors_inner_elastic) :: &
+ num_elem_colors_elastic
+
! local parameters
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: v_tmp
integer,dimension(:),allocatable :: v_tmp_i
@@ -372,6 +384,18 @@
if(num_phase_ispec_poroelastic > 0 ) write(IOUT) phase_ispec_inner_poroelastic
endif
+ ! mesh coloring
+ if( USE_MESH_COLORING_GPU ) then
+ if( ACOUSTIC_SIMULATION ) then
+ write(IOUT) num_colors_outer_acoustic,num_colors_inner_acoustic
+ write(IOUT) num_elem_colors_acoustic
+ endif
+ if( ELASTIC_SIMULATION ) then
+ write(IOUT) num_colors_outer_elastic,num_colors_inner_elastic
+ write(IOUT) num_elem_colors_elastic
+ endif
+ endif
+
close(IOUT)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/check_mesh_quality.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/check_mesh_quality.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/check_mesh_quality.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -54,7 +54,7 @@
logical :: CREATE_VTK_FILES
character(len=256) prname
-
+
! local parameters
integer :: ispec,ispec_min_edge_length,ispec_max_edge_length,ispec_max_skewness, &
ispec_max_skewness_MPI,skewness_max_rank,NSPEC_ALL_SLICES
@@ -124,6 +124,7 @@
ispec_min_edge_length = -1
ispec_max_edge_length = -1
+ ispec_max_skewness = -1
! debug: for vtk output
if( CREATE_VTK_FILES ) then
@@ -146,7 +147,7 @@
if(equiangle_skewness > equiangle_skewness_max) ispec_max_skewness = ispec
if( CREATE_VTK_FILES ) tmp1(ispec) = equiangle_skewness
-
+
! compute minimum and maximum of quality numbers
equiangle_skewness_min = min(equiangle_skewness_min,equiangle_skewness)
edge_aspect_ratio_min = min(edge_aspect_ratio_min,edge_aspect_ratio)
@@ -350,8 +351,8 @@
end if
! debug: for vtk output
- if( CREATE_VTK_FILES ) then
- ! vtk file output
+ if( CREATE_VTK_FILES ) then
+ ! vtk file output
open(66,file=prname(1:len_trim(prname))//'skewness.vtk',status='unknown')
write(66,'(a)') '# vtk DataFile Version 3.1'
write(66,'(a)') 'material model VTK file'
@@ -360,7 +361,7 @@
write(66, '(a,i12,a)') 'POINTS ', nglob, ' float'
do ipoin = 1,nglob
write(66,*) sngl(x(ipoin)),sngl(y(ipoin)),sngl(z(ipoin))
- enddo
+ enddo
write(66,*) ""
! note: indices for vtk start at 0
@@ -384,8 +385,8 @@
write(66,*) tmp1(ispec)
enddo
write(66,*) ""
- close(66)
-
+ close(66)
+
deallocate(tmp1)
endif
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/compute_parameters.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/compute_parameters.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -315,7 +315,7 @@
NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
- !debug
+ !debug
!print*,'nspec minmax:',NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_yMAX
! theoretical number of Gauss-Lobatto points in radial direction
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -397,7 +397,7 @@
if( ier /= 0 ) stop 'error allocating array nodes_coords'
nodes_coords(:,:) = 0.0d0
ibool(:,:,:,:) = 0
-
+
do ispec=1,nspec
ieoff = NGLLCUBE*(ispec-1)
ilocnum = 0
@@ -413,7 +413,7 @@
enddo
enddo
enddo
-
+
! checks ibool range
if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= nglob) then
print*,'error ibool: maximum value ',maxval(ibool(:,:,:,:)) ,'should be ',nglob
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_visual_files.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_visual_files.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/create_visual_files.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -71,7 +71,7 @@
ibool(2,2,2,ispec)
end do
close(64)
-
+
end if
@@ -127,7 +127,7 @@
end if
if( CREATE_VTK_FILES ) then
- ! vtk file output
+ ! vtk file output
open(66,file=prname(1:len_trim(prname))//'mesh.vtk',status='unknown')
write(66,'(a)') '# vtk DataFile Version 3.1'
write(66,'(a)') 'material model VTK file'
@@ -136,7 +136,7 @@
write(66, '(a,i12,a)') 'POINTS ', nglob, ' float'
do ipoin = 1,nglob
write(66,*) sngl(nodes_coords(ipoin,1)),sngl(nodes_coords(ipoin,2)),sngl(nodes_coords(ipoin,3))
- enddo
+ enddo
write(66,*) ""
! note: indices for vtk start at 0
@@ -161,7 +161,7 @@
enddo
write(66,*) ""
close(66)
-
+
endif
call sync_all()
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/meshfem3D.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/meshfem3D.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -405,7 +405,7 @@
write(IMAIN,*) 'error: number of MPI processors actually run on: ',sizeprocs
print*
print*, 'error meshfem3D: number of processors supposed to run on: ',NPROC
- print*, 'error meshfem3D: number of MPI processors actually run on: ',sizeprocs
+ print*, 'error meshfem3D: number of MPI processors actually run on: ',sizeprocs
print*
endif
call exit_MPI(myrank,'wrong number of MPI processes')
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/store_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/store_boundaries.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/meshfem3D/store_boundaries.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -90,21 +90,21 @@
! on boundary: ymax
if(iboun(4,ispec)) then
ispecb4=ispecb4+1
- if( ispecb4 > NSPEC2DMAX_YMIN_YMAX ) stop 'error NSPEC2DMAX_YMIN_YMAX too small'
+ if( ispecb4 > NSPEC2DMAX_YMIN_YMAX ) stop 'error NSPEC2DMAX_YMIN_YMAX too small'
ibelm_ymax(ispecb4)=ispec
endif
! on boundary: bottom
if(iboun(5,ispec)) then
ispecb5=ispecb5+1
- if( ispecb5 > NSPEC2D_BOTTOM ) stop 'error NSPEC2D_BOTTOM too small'
+ if( ispecb5 > NSPEC2D_BOTTOM ) stop 'error NSPEC2D_BOTTOM too small'
ibelm_bottom(ispecb5)=ispec
endif
! on boundary: top
if(iboun(6,ispec)) then
ispecb6=ispecb6+1
- if( ispecb6 > NSPEC2D_TOP ) stop 'error NSPEC2D_TOP too small'
+ if( ispecb6 > NSPEC2D_TOP ) stop 'error NSPEC2D_TOP too small'
ibelm_top(ispecb6)=ispec
endif
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -407,9 +407,11 @@
!real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
logical :: has_vs_zero,has_vp2_zero
+ real(kind=CUSTOM_REAL),dimension(1) :: tmp_val
! debug: for vtk output
real(kind=CUSTOM_REAL),dimension(:),allocatable :: tmp1,tmp2
+
integer:: ier
character(len=256) :: filename,prname
@@ -685,6 +687,7 @@
! returns minimum period
if( myrank == 0 ) min_resolved_period = pmax_glob
+
tmp_val(1) = min_resolved_period
call bcast_all_cr(tmp_val,1)
min_resolved_period = tmp_val(1)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/combine_vol_data.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/combine_vol_data.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -98,7 +98,7 @@
logical :: ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
character(len=256) LOCAL_PATH
integer :: IMODEL
-
+
! checks given arguments
print *
print *,'Recombining ParaView data for slices'
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/constants.h.in 2012-05-13 20:21:06 UTC (rev 20103)
@@ -424,3 +424,4 @@
integer, parameter :: IMODEL_USER_EXTERNAL = 6
integer, parameter :: IMODEL_GLL = 7
integer, parameter :: IMODEL_SALTON_TROUGH = 8
+ integer, parameter :: IMODEL_1D_PREM_PB = 9
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -50,7 +50,7 @@
logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,SUPPRESS_UTM_PROJECTION
character(len=256) LOCAL_PATH,CMTSOLUTION
-
+
! local variables
integer ::ios,icounter,isource,idummy,nproc_eta_old,nproc_xi_old
double precision :: hdur,minval_hdur
@@ -101,7 +101,7 @@
! define the velocity model
call read_value_string(MODEL, 'model.MODEL')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file: MODEL'
-
+
call read_value_logical(OCEANS, 'model.OCEANS')
if(err_occurred() /= 0) return
call read_value_logical(TOPOGRAPHY, 'model.TOPOGRAPHY')
@@ -242,7 +242,7 @@
case( '1d_cascadia')
IMODEL = IMODEL_1D_CASCADIA
- ! user models
+ ! user models
case( 'salton_trough')
IMODEL = IMODEL_SALTON_TROUGH
case( 'tomo' )
@@ -252,15 +252,18 @@
case( 'aniso' )
IMODEL = IMODEL_DEFAULT
ANISOTROPY = .true.
- case default
+ case( '1d_prem_pb' )
+ IMODEL = IMODEL_1D_PREM_PB
+
+ case default
print*
print*,'********** model not recognized: ',trim(MODEL),' **************'
print*,'********** using model: default',' **************'
print*
IMODEL = IMODEL_DEFAULT
end select
-
+
end subroutine read_parameter_file
!
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_topo_bathy_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_topo_bathy_file.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_topo_bathy_file.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -48,7 +48,7 @@
print*,'error opening topography file: ',trim(TOPO_FILE)
stop 'error opening topography file'
endif
-
+
! reads in values
do iy=1,NY_TOPO
do ix=1,NX_TOPO
@@ -74,7 +74,7 @@
include "constants.h"
real(kind=CUSTOM_REAL),intent(in) :: x_target,y_target
-
+
real(kind=CUSTOM_REAL),intent(out) :: target_elevation
integer :: NX_TOPO,NY_TOPO
@@ -87,7 +87,7 @@
double precision :: xval,yval,long,lat
double precision :: long_corner,lat_corner,ratio_xi,ratio_eta
integer :: icornerlong,icornerlat
-
+
! get coordinates of current point
xval = dble(x_target)
yval = dble(y_target)
@@ -127,7 +127,7 @@
itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
end subroutine get_topo_bathy_elevation
-
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -143,10 +143,10 @@
include "constants.h"
real(kind=CUSTOM_REAL),intent(in) :: x_target,y_target
-
+
real(kind=CUSTOM_REAL),intent(out) :: target_elevation
real(kind=CUSTOM_REAL),intent(out) :: target_distmin
-
+
integer :: NSPEC_AB,NGLOB_AB
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -160,7 +160,7 @@
integer, dimension(3,NGLLSQUARE,num_free_surface_faces) :: free_surface_ijk
! local parameters
- real(kind=CUSTOM_REAL),dimension(4) :: elevation_node,dist_node
+ real(kind=CUSTOM_REAL),dimension(4) :: elevation_node,dist_node
real(kind=CUSTOM_REAL) :: distmin,dist
integer :: iface,i,j,ispec,iglob,igll,jgll,kgll
@@ -172,26 +172,26 @@
integer,parameter :: MIDX = (NGLLX+1)/2
integer,parameter :: MIDY = (NGLLY+1)/2
integer,parameter :: MIDZ = (NGLLZ+1)/2
-
+
real(kind=CUSTOM_REAL) :: typical_size
logical :: located_target
-
+
! initialize
target_elevation = 0.0_CUSTOM_REAL
target_distmin = HUGEVAL
-
-
+
+
if(num_free_surface_faces > 0) then
! computes typical size of elements at the surface (uses first element for estimation)
if( USE_DISTANCE_CRITERION ) then
- ispec = free_surface_ispec(1)
+ ispec = free_surface_ispec(1)
typical_size = (xstore(ibool(1,1,1,ispec)) - xstore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 &
+ (ystore(ibool(1,1,1,ispec)) - ystore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2
! use 10 times the distance as a criterion for point detection
typical_size = 10. * typical_size
endif
-
+
! flag to check that we located at least one target element
located_target = .false.
@@ -200,7 +200,7 @@
iselected = 2
jselected = 2
iface_selected = 1
-
+
! loops over all free surface faces
do iface=1,num_free_surface_faces
ispec = free_surface_ispec(iface)
@@ -208,10 +208,10 @@
! exclude elements that are too far from target
if( USE_DISTANCE_CRITERION ) then
iglob = ibool(MIDX,MIDY,MIDZ,ispec)
- dist = (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2
+ dist = (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2
if( dist > typical_size ) cycle
endif
-
+
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
do j = 2,NGLLY - 1
@@ -220,13 +220,13 @@
igll = free_surface_ijk(1,(j-1)*NGLLY+i,iface)
jgll = free_surface_ijk(2,(j-1)*NGLLY+i,iface)
kgll = free_surface_ijk(3,(j-1)*NGLLY+i,iface)
-
+
iglob = ibool(igll,jgll,kgll,ispec)
! distance (squared) to target
dist = ( x_target - xstore(iglob) )**2 + &
( y_target - ystore(iglob) )**2
-
+
! keep this point if it is closer to the receiver
if(dist < distmin) then
distmin = dist
@@ -238,8 +238,8 @@
located_target = .true.
endif
enddo
- enddo
- end do
+ enddo
+ end do
! if we have not located a target element, the point is not in this slice
! therefore use first element only for fictitious iterative search
@@ -248,7 +248,7 @@
jselected = 2
iface_selected = 1
endif
-
+
! weighted mean at current point of topography elevation of the four closest nodes
! set distance to huge initial value
distmin = HUGEVAL
@@ -264,23 +264,23 @@
jgll = free_surface_ijk(2,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
kgll = free_surface_ijk(3,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
iglob = ibool(igll,jgll,kgll,ispec)
-
+
! stores node infos
inode = inode + 1
elevation_node(inode) = zstore(iglob)
dist_node(inode) = sqrt( (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2 )
end do
end do
-
+
! weighted elevation
dist = sum( dist_node(:) )
if(dist < distmin) then
-
+
! sets new minimum distance (of all 4 closest nodes)
distmin = dist
target_distmin = distmin
-
- ! interpolates elevation
+
+ ! interpolates elevation
if( dist > TINYVAL ) then
target_elevation = (dist_node(1)/dist)*elevation_node(1) + &
(dist_node(2)/dist)*elevation_node(2) + &
@@ -290,10 +290,10 @@
stop 'error summed distance to node is zero'
endif
endif
-
- end do
+
+ end do
end do
-
+
end if
end subroutine get_topo_elevation_free
@@ -313,10 +313,10 @@
include "constants.h"
real(kind=CUSTOM_REAL),intent(in) :: x_target,y_target
-
+
real(kind=CUSTOM_REAL),intent(out) :: target_elevation
real(kind=CUSTOM_REAL),intent(out) :: target_distmin
-
+
integer :: NSPEC_AB,NGLOB_AB
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -339,26 +339,26 @@
integer,parameter :: MIDX = (NGLLX+1)/2
integer,parameter :: MIDY = (NGLLY+1)/2
integer,parameter :: MIDZ = (NGLLZ+1)/2
-
+
real(kind=CUSTOM_REAL) :: typical_size
logical :: located_target
-
+
! initialize
target_elevation = 0.0_CUSTOM_REAL
target_distmin = HUGEVAL
-
-
+
+
if(num_free_surface_faces > 0) then
! computes typical size of elements at the surface (uses first element for estimation)
if( USE_DISTANCE_CRITERION ) then
- ispec = free_surface_ispec(1)
+ ispec = free_surface_ispec(1)
typical_size = (xstore(ibool(1,1,1,ispec)) - xstore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2 &
+ (ystore(ibool(1,1,1,ispec)) - ystore(ibool(NGLLX,NGLLY,NGLLZ,ispec)))**2
! use 10 times the distance as a criterion for point detection
typical_size = 10. * typical_size
endif
-
+
! flag to check that we located at least one target element
located_target = .false.
@@ -372,22 +372,22 @@
! excludes elements that are too far from target
if( USE_DISTANCE_CRITERION ) then
iglob = ibool(MIDX,MIDY,MIDZ,ispec)
- dist = (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2
+ dist = (x_target - xstore(iglob))**2 + (y_target - ystore(iglob))**2
if( dist > typical_size ) cycle
endif
-
+
! loop only on points inside the element
do i = 1,NGLLSQUARE
igll = free_surface_ijk(1,i,iface)
jgll = free_surface_ijk(2,i,iface)
kgll = free_surface_ijk(3,i,iface)
-
+
iglob = ibool(igll,jgll,kgll,ispec)
! distance (squared) to target
dist = ( x_target - xstore(iglob) )**2 + &
( y_target - ystore(iglob) )**2
-
+
! keep this point if it is closer to the receiver
if(dist < distmin) then
distmin = dist
@@ -397,8 +397,8 @@
target_distmin = dist
located_target = .true.
endif
- enddo
- end do
+ enddo
+ end do
! if we have not located a target element, the point is not in this slice
! therefore use first element only for fictitious iterative search
@@ -409,10 +409,10 @@
! elevation (given in z - coordinate)
target_elevation = zstore(iglob)
target_distmin = ( x_target - xstore(iglob) )**2 + ( y_target - ystore(iglob) )**2
- located_target = .true.
+ located_target = .true.
endif
-
+
end if
end subroutine get_topo_elevation_free_closest
-
\ No newline at end of file
+
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -38,6 +38,8 @@
!
double precision function wtime()
+
+ implicit none
real :: ct
! note: for simplicity, we take cpu_time which returns the elapsed CPU time in seconds
@@ -45,6 +47,7 @@
call cpu_time(ct)
wtime = ct
+
end function wtime
!
@@ -682,11 +685,10 @@
subroutine send_dp(sendbuf, sendcount, dest, sendtag)
implicit none
- include "constants.h"
integer dest,sendtag
integer sendcount
- real(kind=CUSTOM_REAL),dimension(sendcount):: sendbuf
+ double precision,dimension(sendcount):: sendbuf
stop 'send_dp not implemented for serial code'
@@ -697,13 +699,12 @@
subroutine recv_dp(recvbuf, recvcount, dest, recvtag)
implicit none
- include "constants.h"
integer dest,recvtag
integer recvcount
- real(kind=CUSTOM_REAL),dimension(recvcount):: recvbuf
+ double precision,dimension(recvcount):: recvbuf
- stop 'recv_dp not implemented for parallel code'
+ stop 'recv_dp not implemented for serial code'
end subroutine recv_dp
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/smooth_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/smooth_vol_data.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/smooth_vol_data.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -73,7 +73,7 @@
real(kind=CUSTOM_REAL), dimension(:,:),allocatable :: dummy_2
real(kind=CUSTOM_REAL), dimension(:,:,:),allocatable :: dummy_3
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:),allocatable :: dummy_5
-
+
integer, dimension(:),allocatable :: idummy
logical, dimension(:),allocatable :: ldummy
integer, dimension(:,:,:),allocatable :: idummy_3
@@ -262,7 +262,7 @@
! reads mesh file
!
! needs to get point locations, jacobians and MPI neighbours
-
+
! opens external mesh file
write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'/proc',myrank,'_'//'external_mesh.bin'
open(unit=27,file=trim(prname_lp),&
@@ -273,7 +273,7 @@
call exit_mpi(myrank, 'error reading external mesh file')
endif
- ! gets number of elements and global points for this partition
+ ! gets number of elements and global points for this partition
read(27) NSPEC_AB
read(27) NGLOB_AB
@@ -292,7 +292,7 @@
allocate(dummy(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array dummy and jacobian'
-
+
! needs jacobian
read(27) dummy ! xix
read(27) dummy ! xiy
@@ -306,7 +306,7 @@
read(27) jacobian
! now skips all until MPI section can be read in
-
+
! reads in partiton neighbors
read(27) dummy ! kappastore
read(27) dummy ! mustore
@@ -331,7 +331,7 @@
read(27) dummy_1 ! rmass_acoustic
read(27) dummy ! rhostore
endif
-
+
! elastic
if( ELASTIC_SIMULATION ) then
read(27) dummy_1 ! rmass
@@ -341,28 +341,28 @@
read(27) dummy ! rho_vp
read(27) dummy ! rho_vs
endif
-
+
! poroelastic
if( POROELASTIC_SIMULATION ) then
read(27) dummy_1 ! rmass_solid_poroelastic
read(27) dummy_1 ! rmass_fluid_poroelastic
allocate(dummy_5(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB))
read(27) dummy_5 ! rhoarraystore
- deallocate(dummy_5)
+ deallocate(dummy_5)
allocate(dummy_5(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB))
read(27) dummy_5 ! kappaarraystore
deallocate(dummy_5)
read(27) dummy ! etastore
- read(27) dummy ! tortstore
+ read(27) dummy ! tortstore
allocate(dummy_5(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB))
read(27) dummy_5 ! permstore
deallocate(dummy_5)
read(27) dummy ! phistore
read(27) dummy ! rho_vpI
read(27) dummy ! rho_vpII
- read(27) dummy ! rho_vsI
+ read(27) dummy ! rho_vsI
endif
-
+
deallocate(dummy_1)
deallocate(dummy)
@@ -407,7 +407,7 @@
idummy_3(3,NGLLSQUARE,idummy_a), &
dummy_2(NGLLSQUARE,idummy_a), &
dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
- if( ier /= 0 ) stop 'error allocating array idummy etc.'
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
read(27) idummy ! coupling_ac_el_ispec
read(27) idummy_3 ! coupling_ac_el_ijk
read(27) dummy_2 ! coupling_ac_el_jacobian2Dw
@@ -422,7 +422,7 @@
idummy_3(3,NGLLSQUARE,idummy_a), &
dummy_2(NGLLSQUARE,idummy_a), &
dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
- if( ier /= 0 ) stop 'error allocating array idummy etc.'
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
read(27) idummy ! coupling_ac_po_ispec
read(27) idummy_3 ! coupling_ac_po_ijk
read(27) dummy_2 ! coupling_ac_po_jacobian2Dw
@@ -437,11 +437,11 @@
idummy_3(3,NGLLSQUARE,idummy_a), &
dummy_2(NGLLSQUARE,idummy_a), &
dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
- if( ier /= 0 ) stop 'error allocating array idummy etc.'
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
read(27) idummy ! coupling_el_po_ispec
- read(27) idummy ! coupling_po_el_ispec
+ read(27) idummy ! coupling_po_el_ispec
read(27) idummy_3 ! coupling_el_po_ijk
- read(27) idummy_3 ! coupling_po_el_ijk
+ read(27) idummy_3 ! coupling_po_el_ijk
read(27) dummy_2 ! coupling_el_po_jacobian2Dw
read(27) dummy_3 ! coupling_el_po_normal
deallocate( idummy,idummy_3,dummy_2,dummy_3)
@@ -457,12 +457,12 @@
read(27) my_neighbours_ext_mesh
! no more information is needed from external mesh files
endif
-
+
! we're done reading in mesh arrays
close(27)
! ---------------------
-
+
! for smoothing, we use cell centers to find and locate nearby elements
!
! sets the location of the center of the elements and local points
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in 2012-05-13 20:21:06 UTC (rev 20103)
@@ -41,7 +41,8 @@
# with configure: ./configure --with-cuda CUDA_LIB=.. CUDA_INC=.. MPI_INC=..
# default cuda libraries
- at COND_CUDA_TRUE@CUDA_LIBS = -lcuda -lcudart -lcublas
+# runtime library -lcudart needed, others are optional -lcuda -lcublas
+ at COND_CUDA_TRUE@CUDA_LIBS = -lcudart
@COND_CUDA_FALSE at CUDA_LIBS =
CUDA_LIB_LOCATION = @CUDA_LIB@
@@ -214,19 +215,18 @@
$O/model_update.o \
$O/specfem3D_par.o \
$O/initialize_simulation.o \
- $O/read_parameter_file.o \
- $O/read_value_parameters.o \
- $O/get_value_parameters.o \
- $O/param_reader.o \
- $O/exit_mpi.o \
- $O/parallel.o \
+ $O/read_parameter_file.shared.o \
+ $O/read_value_parameters.shared.o \
+ $O/get_value_parameters.shared.o \
+ $O/param_reader.cc.o \
+ $O/exit_mpi.shared.o \
$O/read_mesh_databases.o \
- $O/create_name_database.o \
- $O/check_mesh_resolution.o \
- $O/gll_library.o \
- $O/get_attenuation_model.o \
+ $O/create_name_database.shared.o \
+ $O/check_mesh_resolution.shared.o \
+ $O/gll_library.shared.o \
+ $O/get_attenuation_model.shared.o \
$O/save_external_bin_m_up.o \
- $O/write_VTK_data.o \
+ $O/write_VTK_data.shared.o \
$(EMPTY_MACRO)
# objects toggled between the parallel and serial version
@@ -297,8 +297,8 @@
xsmooth_vol_data: $O/smooth_vol_data.o $O/read_parameter_file.shared.o $O/read_value_parameters.shared.o $O/get_value_parameters.shared.o $O/param_reader.cc.o $O/gll_library.shared.o $O/exit_mpi.shared.o $(COND_MPI_OBJECTS)
${FCLINK} -o ${E}/xsmooth_vol_data $O/smooth_vol_data.o $O/read_parameter_file.shared.o $O/read_value_parameters.shared.o $O/get_value_parameters.shared.o $O/param_reader.cc.o $O/gll_library.shared.o $O/exit_mpi.shared.o $(COND_MPI_OBJECTS) $(MPILIBS)
-xmodel_update: $(MODEL_UPD_OBJECTS)
- ${FCLINK} -o ${E}/xmodel_update $(MODEL_UPD_OBJECTS) $(MPILIBS)
+xmodel_update: $(MODEL_UPD_OBJECTS) $(COND_MPI_OBJECTS)
+ ${FCLINK} -o ${E}/xmodel_update $(MODEL_UPD_OBJECTS) $(COND_MPI_OBJECTS) $(MPILIBS)
clean:
rm -f $O/* *.o *.gnu *.mod $(OUTPUT)/timestamp* $(OUTPUT)/starttime*txt work.pc* \
@@ -335,72 +335,21 @@
$O/%.shared.o: ${SHARED}%.F90 $(SHARED)constants.h
${FCCOMPILE_NO_CHECK} -c -o $@ $<
+
###
### OpenMP compilation
###
$O/%.openmp.o: %.f90 $(SHARED)constants.h
${FCCOMPILE_NO_CHECK} -c -o $@ $<
-$O/compute_forces_poroelastic.o: $(SHARED)constants.h compute_forces_poroelastic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_poroelastic.o compute_forces_poroelastic.f90
-$O/compute_forces_solid.o: $(SHARED)constants.h compute_forces_solid.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_solid.o compute_forces_solid.f90
+###
+### CUDA compilation
+###
+$O/%.cuda.o: ${CUDAD}%.cu ../../config.h $(CUDAD)mesh_constants_cuda.h $(CUDAD)prepare_constants_cuda.h
+ $(NVCC) -c $< -o $@ $(NVCC_FLAGS)
-$O/compute_forces_fluid.o: $(SHARED)constants.h compute_forces_fluid.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_fluid.o compute_forces_fluid.f90
-
-$O/compute_forces_acoustic.o: $(SHARED)constants.h compute_forces_acoustic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_acoustic.o compute_forces_acoustic.f90
-
-$O/compute_forces_acoustic_pot.o: $(SHARED)constants.h compute_forces_acoustic_pot.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_acoustic_pot.o compute_forces_acoustic_pot.f90
-
-$O/compute_forces_acoustic_PML.o: $(SHARED)constants.h compute_forces_acoustic_PML.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_acoustic_PML.o compute_forces_acoustic_PML.f90
-
-$O/compute_add_sources_acoustic.o: $(SHARED)constants.h compute_add_sources_acoustic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_add_sources_acoustic.o compute_add_sources_acoustic.f90
-
-$O/compute_add_sources_elastic.o: $(SHARED)constants.h compute_add_sources_elastic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_add_sources_elastic.o compute_add_sources_elastic.f90
-
-$O/compute_add_sources_poroelastic.o: $(SHARED)constants.h compute_add_sources_poroelastic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_add_sources_poroelastic.o compute_add_sources_poroelastic.f90
-
-$O/get_poroelastic_velocities.o: $(SHARED)constants.h get_poroelastic_velocities.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/get_poroelastic_velocities.o get_poroelastic_velocities.f90
-
-$O/compute_coupling_acoustic_el.o: $(SHARED)constants.h compute_coupling_acoustic_el.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_acoustic_el.o compute_coupling_acoustic_el.f90
-
-$O/compute_coupling_elastic_ac.o: $(SHARED)constants.h compute_coupling_elastic_ac.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_elastic_ac.o compute_coupling_elastic_ac.f90
-
-$O/compute_coupling_elastic_po.o: $(SHARED)constants.h compute_coupling_elastic_po.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_elastic_po.o compute_coupling_elastic_po.f90
-
-$O/compute_coupling_acoustic_po.o: $(SHARED)constants.h compute_coupling_acoustic_po.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_acoustic_po.o compute_coupling_acoustic_po.f90
-
-$O/compute_coupling_poroelastic_ac.o: $(SHARED)constants.h compute_coupling_poroelastic_ac.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_poroelastic_ac.o compute_coupling_poroelastic_ac.f90
-
-$O/compute_coupling_poroelastic_el.o: $(SHARED)constants.h compute_coupling_poroelastic_el.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_poroelastic_el.o compute_coupling_poroelastic_el.f90
-
-$O/compute_stacey_acoustic.o: $(SHARED)constants.h compute_stacey_acoustic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_stacey_acoustic.o compute_stacey_acoustic.f90
-
-$O/compute_stacey_elastic.o: $(SHARED)constants.h compute_stacey_elastic.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_stacey_elastic.o compute_stacey_elastic.f90
-
-$O/compute_gradient.o: $(SHARED)constants.h compute_gradient.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_gradient.o compute_gradient.f90
-
-$O/compute_interpolated_dva.o: $(SHARED)constants.h compute_interpolated_dva.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_interpolated_dva.o compute_interpolated_dva.f90
-
+###
### C compilation
###
force_ftz.o: ${SHARED}/force_ftz.c ../../config.h
@@ -413,55 +362,6 @@
${CC} -c $(CFLAGS) $(MPI_INC) -o $@ ${CUDAD}/$< -I../../
-$O/detect_mesh_surfaces.o: $(SHARED)constants.h detect_mesh_surfaces.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/detect_mesh_surfaces.o detect_mesh_surfaces.f90
-
-$O/setup_movie_meshes.o: $(SHARED)constants.h setup_movie_meshes.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/setup_movie_meshes.o setup_movie_meshes.f90
-
-$O/read_topography_bathymetry.o: $(SHARED)constants.h read_topography_bathymetry.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/read_topography_bathymetry.o read_topography_bathymetry.f90
-
-$O/setup_sources_receivers.o: $(SHARED)constants.h setup_sources_receivers.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/setup_sources_receivers.o setup_sources_receivers.f90
-
-$O/prepare_timerun.o: $(SHARED)constants.h prepare_timerun.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/prepare_timerun.o prepare_timerun.f90
-
-$O/iterate_time.o: $(SHARED)constants.h iterate_time.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/iterate_time.o iterate_time.f90
-
-$O/finalize_simulation.o: $(SHARED)constants.h finalize_simulation.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/finalize_simulation.o finalize_simulation.f90
-
-$O/assemble_MPI_vector.o: $(SHARED)constants.h assemble_MPI_vector.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_vector.o assemble_MPI_vector.f90
-
-$O/assemble_MPI_scalar.o: $(SHARED)constants.h ${SHARED}/assemble_MPI_scalar.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_scalar.o ${SHARED}/assemble_MPI_scalar.f90
-
-$O/save_adjoint_kernels.o: $(SHARED)constants.h save_adjoint_kernels.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/save_adjoint_kernels.o save_adjoint_kernels.f90
-
-$O/write_movie_output.o: $(SHARED)constants.h write_movie_output.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/write_movie_output.o write_movie_output.f90
-
-$O/create_color_image.o: $(SHARED)constants.h create_color_image.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/create_color_image.o create_color_image.f90
-
-$O/write_seismograms.o: $(SHARED)constants.h write_seismograms.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/write_seismograms.o write_seismograms.f90
-
-$O/write_output_ASCII.o: $(SHARED)constants.h write_output_ASCII.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/write_output_ASCII.o write_output_ASCII.f90
-
-$O/write_output_SU.o: $(SHARED)constants.h write_output_SU.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/write_output_SU.o write_output_SU.f90
-
-$O/noise_tomography.o: $(SHARED)constants.h noise_tomography.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/noise_tomography.o noise_tomography.f90
-
-
###
### MPI compilation without optimization
###
@@ -471,138 +371,14 @@
$O/serial.o: $(SHARED)constants.h ${SHARED}/serial.f90
${FCCOMPILE_CHECK} -c -o $O/serial.o ${SHARED}/serial.f90
+
+$O/smooth_vol_data.o: $(SHARED)constants.h ${SHARED}/smooth_vol_data.f90
+ ${MPIFCCOMPILE_NO_CHECK} -c -o $O/smooth_vol_data.o ${SHARED}/smooth_vol_data.f90
-$O/program_specfem3D.o: program_specfem3D.f90
- ${FCCOMPILE_CHECK} -c -o $O/program_specfem3D.o program_specfem3D.f90
-$O/locate_source.o: $(SHARED)constants.h locate_source.f90
- ${FCCOMPILE_CHECK} -c -o $O/locate_source.o locate_source.f90
-$O/locate_receivers.o: $(SHARED)constants.h locate_receivers.f90
- ${FCCOMPILE_CHECK} -c -o $O/locate_receivers.o locate_receivers.f90
-$O/exit_mpi.o: $(SHARED)constants.h ${SHARED}/exit_mpi.f90
- ${FCCOMPILE_CHECK} -c -o $O/exit_mpi.o ${SHARED}/exit_mpi.f90
-
-$O/convolve_source_timefunction.o: ${SHARED}/convolve_source_timefunction.f90
- ${FCCOMPILE_CHECK} -c -o $O/convolve_source_timefunction.o ${SHARED}/convolve_source_timefunction.f90
-
-$O/save_header_file.o: $(SHARED)constants.h ${SHARED}/save_header_file.f90
- ${FCCOMPILE_CHECK} -c -o $O/save_header_file.o ${SHARED}/save_header_file.f90
-
-$O/read_parameter_file.o: $(SHARED)constants.h ${SHARED}/read_parameter_file.f90
- ${FCCOMPILE_CHECK} -c -o $O/read_parameter_file.o ${SHARED}/read_parameter_file.f90
-
-$O/read_value_parameters.o: $(SHARED)constants.h ${SHARED}/read_value_parameters.f90
- ${FCCOMPILE_CHECK} -c -o $O/read_value_parameters.o ${SHARED}/read_value_parameters.f90
-
-$O/get_value_parameters.o: $(SHARED)constants.h ${SHARED}/get_value_parameters.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_value_parameters.o ${SHARED}/get_value_parameters.f90
-
-$O/utm_geo.o: $(SHARED)constants.h ${SHARED}/utm_geo.f90
- ${FCCOMPILE_CHECK} -c -o $O/utm_geo.o ${SHARED}/utm_geo.f90
-
-$O/check_mesh_resolution.o: $(SHARED)constants.h ${SHARED}/check_mesh_resolution.f90
- ${FCCOMPILE_CHECK} -c -o $O/check_mesh_resolution.o ${SHARED}/check_mesh_resolution.f90
-
-$O/detect_surface.o: $(SHARED)constants.h ${SHARED}/detect_surface.f90
- ${FCCOMPILE_CHECK} -c -o $O/detect_surface.o ${SHARED}/detect_surface.f90
-
-$O/gll_library.o: $(SHARED)constants.h ${SHARED}/gll_library.f90
- ${FCCOMPILE_CHECK} -c -o $O/gll_library.o ${SHARED}/gll_library.f90
-
-$O/get_jacobian_boundaries.o: $(SHARED)constants.h ${SHARED}/get_jacobian_boundaries.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_jacobian_boundaries.o ${SHARED}/get_jacobian_boundaries.f90
-
-$O/get_flags_boundaries.o: $(SHARED)constants.h ${SHARED}/get_flags_boundaries.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_flags_boundaries.o ${SHARED}/get_flags_boundaries.f90
-
-$O/get_cmt.o: $(SHARED)constants.h ${SHARED}/get_cmt.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_cmt.o ${SHARED}/get_cmt.f90
-
-$O/create_movie_shakemap_AVS_DX_GMT.o: $(SHARED)constants.h ${SHARED}/create_movie_shakemap_AVS_DX_GMT.f90 $(OUTPUT)/surface_from_mesher.h
- ${FCCOMPILE_CHECK} -c -o $O/create_movie_shakemap_AVS_DX_GMT.o ${SHARED}/create_movie_shakemap_AVS_DX_GMT.f90 -I$(OUTPUT)
-
-$O/get_element_face.o: $(SHARED)constants.h ${SHARED}/get_element_face.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_element_face.o ${SHARED}/get_element_face.f90
-
-$O/write_VTK_data.o: $(SHARED)constants.h ${SHARED}/write_VTK_data.f90
- ${FCCOMPILE_CHECK} -c -o $O/write_VTK_data.o ${SHARED}/write_VTK_data.f90
-
-$O/get_shape3D.o: $(SHARED)constants.h ${SHARED}/get_shape3D.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_shape3D.o ${SHARED}/get_shape3D.f90
-
-$O/get_shape2D.o: $(SHARED)constants.h ${SHARED}/get_shape2D.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_shape2D.o ${SHARED}/get_shape2D.f90
-
-$O/hex_nodes.o: $(SHARED)constants.h ${SHARED}/hex_nodes.f90
- ${FCCOMPILE_CHECK} -c -o $O/hex_nodes.o ${SHARED}/hex_nodes.f90
-
-$O/netlib_specfun_erf.o: ${SHARED}/netlib_specfun_erf.f90
- ${FCCOMPILE_CHECK} -c -o $O/netlib_specfun_erf.o ${SHARED}/netlib_specfun_erf.f90
-
-$O/sort_array_coordinates.o: $(SHARED)constants.h ${SHARED}/sort_array_coordinates.f90
- ${FCCOMPILE_CHECK} -c -o $O/sort_array_coordinates.o ${SHARED}/sort_array_coordinates.f90
-
-$O/comp_source_time_function.o: $(SHARED)constants.h comp_source_time_function.f90
- ${FCCOMPILE_CHECK} -c -o $O/comp_source_time_function.o comp_source_time_function.f90
-
-$O/read_topo_bathy_file.o: $(SHARED)constants.h ${SHARED}/read_topo_bathy_file.f90
- ${FCCOMPILE_CHECK} -c -o $O/read_topo_bathy_file.o ${SHARED}/read_topo_bathy_file.f90
-
-$O/lagrange_poly.o: $(SHARED)constants.h ${SHARED}/lagrange_poly.f90
- ${FCCOMPILE_CHECK} -c -o $O/lagrange_poly.o ${SHARED}/lagrange_poly.f90
-
-$O/recompute_jacobian.o: $(SHARED)constants.h ${SHARED}/recompute_jacobian.f90
- ${FCCOMPILE_CHECK} -c -o $O/recompute_jacobian.o ${SHARED}/recompute_jacobian.f90
-
-$O/create_name_database.o: $(SHARED)constants.h ${SHARED}/create_name_database.f90
- ${FCCOMPILE_CHECK} -c -o $O/create_name_database.o ${SHARED}/create_name_database.f90
-
-$O/create_serial_name_database.o: $(SHARED)constants.h ${SHARED}/create_serial_name_database.f90
- ${FCCOMPILE_CHECK} -c -o $O/create_serial_name_database.o ${SHARED}/create_serial_name_database.f90
-
-$O/define_derivation_matrices.o: $(SHARED)constants.h ${SHARED}/define_derivation_matrices.f90
- ${FCCOMPILE_CHECK} -c -o $O/define_derivation_matrices.o ${SHARED}/define_derivation_matrices.f90
-
-$O/compute_adj_source_frechet.o: $(SHARED)constants.h compute_adj_source_frechet.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_adj_source_frechet.o compute_adj_source_frechet.f90
-
-$O/compute_arrays_source.o: $(SHARED)constants.h compute_arrays_source.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_arrays_source.o compute_arrays_source.f90
-
-$O/multiply_arrays_source.o: ${SHARED}/constants.h ${SHARED}/multiply_arrays_source.f90
- ${FCCOMPILE_CHECK} -c -o $O/multiply_arrays_source.o ${SHARED}/multiply_arrays_source.f90
-
-$O/get_attenuation_model.o: $(SHARED)constants.h ${SHARED}/get_attenuation_model.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_attenuation_model.o ${SHARED}/get_attenuation_model.f90
-
-$O/compute_boundary_kernel.o: $(SHARED)constants.h compute_boundary_kernel.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o compute_boundary_kernel.f90
-
-$O/compute_kernels.o: $(SHARED)constants.h compute_kernels.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o compute_kernels.f90
-
-$O/combine_vol_data.o: $(SHARED)constants.h ${SHARED}/combine_vol_data.f90
- ${FCCOMPILE_CHECK} -c -o $O/combine_vol_data.o ${SHARED}/combine_vol_data.f90
-
-$O/combine_surf_data.o: $(SHARED)constants.h ${SHARED}/combine_surf_data.f90
- ${FCCOMPILE_CHECK} -c -o $O/combine_surf_data.o ${SHARED}/combine_surf_data.f90
-
-$O/prepare_assemble_MPI.o: $(SHARED)constants.h ${SHARED}/prepare_assemble_MPI.f90
- ${FCCOMPILE_CHECK} -c -o $O/prepare_assemble_MPI.o ${SHARED}/prepare_assemble_MPI.f90
-
-$O/PML_init.o: $(SHARED)constants.h PML_init.f90
- ${FCCOMPILE_CHECK} -c -o $O/PML_init.o PML_init.f90
-
##
-## smoothing
-##
-
-$O/smooth_vol_data.o: $(SHARED)constants.h ${SHARED}/smooth_vol_data.f90
- ${MPIFCCOMPILE_NO_CHECK} -c -o $O/smooth_vol_data.o ${SHARED}/smooth_vol_data.f90
-
-##
## model update
##
@@ -612,13 +388,3 @@
$O/save_external_bin_m_up.o: ${SHARED}/constants.h save_external_bin_m_up.f90
${FCCOMPILE_CHECK} -c -o $O/save_external_bin_m_up.o save_external_bin_m_up.f90
-###
-### C files below
-###
-
-$O/param_reader.o: ${SHARED}/param_reader.c
- ${CC} -c $(CFLAGS) -o $O/param_reader.o ${SHARED}/param_reader.c -I../../
-
-$O/write_c_binary.o: ${SHARED}/write_c_binary.c
- ${CC} -c $(CFLAGS) -o $O/write_c_binary.o ${SHARED}/write_c_binary.c -I../../
-
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/assemble_MPI_vector.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/assemble_MPI_vector.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/assemble_MPI_vector.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -200,62 +200,6 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine assemble_MPI_vector_send_cuda(NPROC, &
- buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh, &
- my_neighbours_ext_mesh, &
- request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
-
- ! sends data
- ! note: array to assemble already filled into buffer_send_vector_ext_mesh array
-
- implicit none
-
- include "constants.h"
-
- integer :: NPROC
-
- integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
-
- real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
- buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
-
- integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
- integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
-
- integer iinterface
-
- ! note: preparation of the contribution between partitions using MPI
- ! already done in transfer_boun_accel routine
-
- ! send only if more than one partition
- if(NPROC > 1) then
-
- ! send messages
- do iinterface = 1, num_interfaces_ext_mesh
- call isend_cr(buffer_send_vector_ext_mesh(1,1,iinterface), &
- NDIM*nibool_interfaces_ext_mesh(iinterface), &
- my_neighbours_ext_mesh(iinterface), &
- itag, &
- request_send_vector_ext_mesh(iinterface) &
- )
- call irecv_cr(buffer_recv_vector_ext_mesh(1,1,iinterface), &
- NDIM*nibool_interfaces_ext_mesh(iinterface), &
- my_neighbours_ext_mesh(iinterface), &
- itag, &
- request_recv_vector_ext_mesh(iinterface) &
- )
- enddo
-
- endif
-
- end subroutine assemble_MPI_vector_send_cuda
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
subroutine assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_AB,array_val, &
buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
@@ -312,74 +256,7 @@
end subroutine assemble_MPI_vector_ext_mesh_w
-!
-!-------------------------------------------------------------------------------------------------
-!
- subroutine assemble_MPI_vector_write_cuda(NPROC,NGLOB_AB,array_val, Mesh_pointer, &
- buffer_recv_vector_ext_mesh, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
- request_send_vector_ext_mesh,request_recv_vector_ext_mesh, &
- FORWARD_OR_ADJOINT )
-
-! waits for data to receive and assembles
-
- implicit none
-
- include "constants.h"
-
- integer :: NPROC
- integer :: NGLOB_AB
- integer(kind=8) :: Mesh_pointer
-! array to assemble
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val
-
- integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
-
- real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
- buffer_recv_vector_ext_mesh
-
- integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
- integer, dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
- integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
-
- integer iinterface ! ipoin
- integer FORWARD_OR_ADJOINT
-
-! here we have to assemble all the contributions between partitions using MPI
-
-! assemble only if more than one partition
- if(NPROC > 1) then
-
-! wait for communications completion (recv)
- do iinterface = 1, num_interfaces_ext_mesh
- call wait_req(request_recv_vector_ext_mesh(iinterface))
- enddo
-
-! adding contributions of neighbours
- call transfer_asmbl_accel_to_device(Mesh_pointer, array_val, buffer_recv_vector_ext_mesh, &
- num_interfaces_ext_mesh, max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,&
- ibool_interfaces_ext_mesh,FORWARD_OR_ADJOINT)
-
- ! This step is done via previous function transfer_and_assemble...
- ! do iinterface = 1, num_interfaces_ext_mesh
- ! do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
- ! array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
- ! array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
- ! enddo
- ! enddo
-
-! wait for communications completion (send)
- do iinterface = 1, num_interfaces_ext_mesh
- call wait_req(request_send_vector_ext_mesh(iinterface))
- enddo
-
- endif
-
- end subroutine assemble_MPI_vector_ext_mesh_w
-
!
!--------------------------------------------------------------------------------------------------
!
@@ -529,3 +406,263 @@
endif
end subroutine assemble_MPI_vector_poro_w
+
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+! with cuda functions...
+
+ subroutine assemble_MPI_vector_send_cuda(NPROC, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
+
+ ! sends data
+ ! note: array to assemble already filled into buffer_send_vector_ext_mesh array
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NPROC
+
+ integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
+
+ integer iinterface
+
+ ! note: preparation of the contribution between partitions using MPI
+ ! already done in transfer_boun_accel routine
+
+ ! send only if more than one partition
+ if(NPROC > 1) then
+
+ ! send messages
+ do iinterface = 1, num_interfaces_ext_mesh
+ call isend_cr(buffer_send_vector_ext_mesh(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_vector_ext_mesh(iinterface) &
+ )
+ call irecv_cr(buffer_recv_vector_ext_mesh(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_vector_ext_mesh(iinterface) &
+ )
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_send_cuda
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine assemble_MPI_vector_write_cuda(NPROC,NGLOB_AB,array_val, Mesh_pointer, &
+ buffer_recv_vector_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh, &
+ FORWARD_OR_ADJOINT )
+
+! waits for data to receive and assembles
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+ integer(kind=8) :: Mesh_pointer
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val
+
+ integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
+ buffer_recv_vector_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
+
+ integer iinterface ! ipoin
+ integer FORWARD_OR_ADJOINT
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! wait for communications completion (recv)
+ do iinterface = 1, num_interfaces_ext_mesh
+ call wait_req(request_recv_vector_ext_mesh(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ call transfer_asmbl_accel_to_device(Mesh_pointer, array_val, buffer_recv_vector_ext_mesh, &
+ num_interfaces_ext_mesh, max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,&
+ ibool_interfaces_ext_mesh,FORWARD_OR_ADJOINT)
+
+ ! This step is done via previous function transfer_and_assemble...
+ ! do iinterface = 1, num_interfaces_ext_mesh
+ ! do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ ! array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ ! array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
+ ! enddo
+ ! enddo
+
+! wait for communications completion (send)
+ do iinterface = 1, num_interfaces_ext_mesh
+ call wait_req(request_send_vector_ext_mesh(iinterface))
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_write_cuda
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine assemble_MPI_scalar_send_cuda(NPROC, &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
+
+! non-blocking MPI send
+
+ ! sends data
+ ! note: assembling data already filled into buffer_send_scalar_ext_mesh array
+ implicit none
+
+ include "constants.h"
+
+ integer :: NPROC
+ integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+
+ real(kind=CUSTOM_REAL), dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
+
+ integer iinterface
+
+! sends only if more than one partition
+ if(NPROC > 1) then
+
+ ! note: partition border copy into the buffer has already been done
+ ! by routine transfer_boun_pot_from_device()
+
+ ! send messages
+ do iinterface = 1, num_interfaces_ext_mesh
+ ! non-blocking synchronous send request
+ call isend_cr(buffer_send_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_scalar_ext_mesh(iinterface) &
+ )
+ ! receive request
+ call irecv_cr(buffer_recv_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_scalar_ext_mesh(iinterface) &
+ )
+
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_scalar_send_cuda
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine assemble_MPI_scalar_write_cuda(NPROC,NGLOB_AB,array_val, &
+ Mesh_pointer, &
+ buffer_recv_scalar_ext_mesh,num_interfaces_ext_mesh, &
+ max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh, &
+ FORWARD_OR_ADJOINT)
+
+! waits for send/receiver to be completed and assembles contributions
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+ integer(kind=8) :: Mesh_pointer
+
+ integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: array_val
+
+
+ real(kind=CUSTOM_REAL), dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
+ buffer_recv_scalar_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
+
+ integer FORWARD_OR_ADJOINT
+
+ integer iinterface ! ipoin
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+ ! wait for communications completion (recv)
+ do iinterface = 1, num_interfaces_ext_mesh
+ call wait_req(request_recv_scalar_ext_mesh(iinterface))
+ enddo
+
+ ! adding contributions of neighbours
+ call transfer_asmbl_pot_to_device(Mesh_pointer, array_val, buffer_recv_scalar_ext_mesh, &
+ num_interfaces_ext_mesh, max_nibool_interfaces_ext_mesh, nibool_interfaces_ext_mesh,&
+ ibool_interfaces_ext_mesh,FORWARD_OR_ADJOINT)
+
+ ! note: adding contributions of neighbours has been done just above for cuda
+ !do iinterface = 1, num_interfaces_ext_mesh
+ ! do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ ! array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ ! array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) &
+ ! + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
+ ! enddo
+ !enddo
+
+ ! wait for communications completion (send)
+ do iinterface = 1, num_interfaces_ext_mesh
+ call wait_req(request_send_scalar_ext_mesh(iinterface))
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_scalar_write_cuda
+
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -41,7 +41,7 @@
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
station_name,network_name,adj_source_file,nrec_local,number_receiver_global, &
- pm1_source_encoding
+ pm1_source_encoding,nsources_local
implicit none
include "constants.h"
@@ -173,18 +173,15 @@
! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
!endif
- ! source encoding
- stf_used = stf_used * pm1_source_encoding(isource)
-
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
- ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
-
! we use nu_source(:,3) here because we want a source normal to the surface.
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
stf_used = FACTOR_FORCE_SOURCE * &
comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+
+ ! source encoding
+ stf_used = stf_used * pm1_source_encoding(isource)
+
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure:
@@ -196,15 +193,6 @@
nint(eta_source(isource)), &
nint(gamma_source(isource)),ispec)
- ! quasi-heaviside
- !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
-
- ! source encoding
- stf = stf * pm1_source_encoding(isource)
-
- ! distinguishes between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
else
if( USE_RICKER_IPATI ) then
@@ -219,6 +207,9 @@
! quasi-heaviside
!stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! source encoding
+ stf = stf * pm1_source_encoding(isource)
+
! distinguishes between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
@@ -247,7 +238,7 @@
stf_used_total = stf_used_total + stf_used
- endif ! ispec_is_elastic
+ endif ! ispec_is_acoustic
endif ! ispec_is_inner
endif ! myrank
@@ -327,78 +318,42 @@
endif
enddo
else
- !!! read SU adjoint sources
- ! range of the block we need to read
- it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
- it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
- write(procname,"(i4)") myrank
- open(unit=IIN_SU1, &
- file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
- access='direct',recl=240+4*(NSTEP))
- do irec_local = 1,nrec_local
- irec = number_receiver_global(irec_local)
- read(IIN_SU1,rec=irec_local) r4head, adj_temp
- adj_src(:,1)=adj_temp(it_start:it_end)
- adj_src(:,2)=0.0 !TRIVIAL
- adj_src(:,3)=0.0 !TRIVIAL
- ! lagrange interpolators for receiver location
- call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
- call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
- call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
- ! interpolates adjoint source onto GLL points within this element
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
- enddo
- enddo
- enddo
- do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
- adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
- enddo
-
- endif
- enddo
- else
- !!! read SU adjoint sources
- ! range of the block we need to read
- it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
- it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
- write(procname,"(i4)") myrank
- open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
+ !!! read SU adjoint sources
+ ! range of the block we need to read
+ it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
+ it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
+ write(procname,"(i4)") myrank
+ open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
status='old',access='direct',recl=240+4*(NSTEP),iostat = ier)
- if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
+ if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj does not exit')
- do irec_local = 1,nrec_local
- irec = number_receiver_global(irec_local)
- read(IIN_SU1,rec=irec_local) r4head, adj_temp
- adj_src(:,1)=adj_temp(it_start:it_end)
- adj_src(:,2)=0.0 !TRIVIAL
- adj_src(:,3)=0.0 !TRIVIAL
- ! lagrange interpolators for receiver location
- call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
- call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
- call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
- ! interpolates adjoint source onto GLL points within this element
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
- enddo
- enddo
- enddo
- do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
- adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
- enddo
- enddo
- close(IIN_SU1)
- endif !if (.not. SU_FORMAT)
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ read(IIN_SU1,rec=irec_local) r4head, adj_temp
+ adj_src(:,1)=adj_temp(it_start:it_end)
+ adj_src(:,2)=0.0 !TRIVIAL
+ adj_src(:,3)=0.0 !TRIVIAL
+ ! lagrange interpolators for receiver location
+ call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
+ call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
+ ! interpolates adjoint source onto GLL points within this element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
+ enddo
+ enddo
+ enddo
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+ enddo
+ close(IIN_SU1)
+ endif !if (.not. SU_FORMAT)
- deallocate(adj_sourcearray)
-
deallocate(adj_sourcearray)
-
endif ! if(ibool_read_adj_arrays)
if( it < NSTEP ) then
@@ -506,13 +461,8 @@
nint(gamma_source(isource)), &
ispec)
- ! source encoding
- stf_used = stf_used * pm1_source_encoding(isource)
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
- ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
-
!if (it == 1 .and. myrank == 0) then
! write(IMAIN,*) 'using a source of dominant frequency ',f0
! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
@@ -527,15 +477,20 @@
stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr( &
dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
- ! quasi-heaviside
- !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! source encoding
+ stf_used = stf_used * pm1_source_encoding(isource)
- ! source encoding
- stf = stf * pm1_source_encoding(isource)
+ ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
+ ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+ ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
- ! distinguishes between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
+ ! acoustic source for pressure gets divided by kappa
+ ! source contribution
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ - stf_used / kappastore(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)),ispec)
+
else
if( USE_RICKER_IPATI ) then
@@ -547,6 +502,12 @@
dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
endif
+ ! quasi-heaviside
+ !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+
+ ! source encoding
+ stf = stf * pm1_source_encoding(isource)
+
! distinguishes between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -269,55 +269,77 @@
allocate(adj_sourcearray(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating array adj_sourcearray'
- endif
- enddo
- else
- !!! read SU adjoint sources
- ! range of the block we need to read
- it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
- it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
- write(procname,"(i4)") myrank
- ! read adjoint sources
- open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
+ if (.not. SU_FORMAT) then
+ !!! read ascii adjoint sources
+ irec_local = 0
+ do irec = 1, nrec
+ ! compute source arrays
+ if (myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
+ ! reads in **sta**.**net**.**LH**.adj files
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+ call compute_arrays_adjoint_source(myrank,adj_source_file, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ adj_sourcearray, xigll,yigll,zigll, &
+ it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+ endif
+ enddo
+ else
+ !!! read SU adjoint sources
+ ! range of the block we need to read
+ it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
+ it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
+ write(procname,"(i4)") myrank
+ ! read adjoint sources
+ open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
status='old',access='direct',recl=240+4*NSTEP,iostat = ier)
- if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
+ if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj does not exit')
- open(unit=IIN_SU2, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dy_SU.adj', &
+ open(unit=IIN_SU2, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dy_SU.adj', &
status='old',access='direct',recl=240+4*NSTEP,iostat = ier)
- if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
+ if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
//'../SEM/'//trim(adjustl(procname))//'_dy_SU.adj does not exit')
- open(unit=IIN_SU3, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dz_SU.adj', &
+ open(unit=IIN_SU3, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dz_SU.adj', &
status='old',access='direct',recl=240+4*NSTEP,iostat = ier)
- if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
+ if( ier /= 0 ) call exit_MPI(myrank,'file '//trim(adjustl(OUTPUT_FILES_PATH)) &
//'../SEM/'//trim(adjustl(procname))//'_dz_SU.adj does not exit')
- do irec_local = 1,nrec_local
- irec = number_receiver_global(irec_local)
- read(IIN_SU1,rec=irec_local) r4head, adj_temp
- adj_src(:,1)=adj_temp(it_start:it_end)
- read(IIN_SU2,rec=irec_local) r4head, adj_temp
- adj_src(:,2)=adj_temp(it_start:it_end)
- read(IIN_SU3,rec=irec_local) r4head, adj_temp
- adj_src(:,3)=adj_temp(it_start:it_end)
- ! lagrange interpolators for receiver location
- call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
- call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
- call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
- ! interpolates adjoint source onto GLL points within this element
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
- enddo
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ read(IIN_SU1,rec=irec_local) r4head, adj_temp
+ adj_src(:,1)=adj_temp(it_start:it_end)
+ read(IIN_SU2,rec=irec_local) r4head, adj_temp
+ adj_src(:,2)=adj_temp(it_start:it_end)
+ read(IIN_SU3,rec=irec_local) r4head, adj_temp
+ adj_src(:,3)=adj_temp(it_start:it_end)
+ ! lagrange interpolators for receiver location
+ call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
+ call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
+ ! interpolates adjoint source onto GLL points within this element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
+ enddo
+ enddo
+ enddo
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+ enddo
+ close(IIN_SU1)
+ close(IIN_SU2)
+ close(IIN_SU3)
+ endif !if(.not. SU_FORMAT)
- deallocate(adj_sourcearray)
-
- endif ! if(ibool_read_adj_arrays)
-
deallocate(adj_sourcearray)
-
endif ! if(ibool_read_adj_arrays)
+
if( it < NSTEP ) then
if(.NOT. GPU_MODE) then
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_poroelastic.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_poroelastic.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -240,7 +240,7 @@
i = adj_sourcearrays(1,1,1,1,1,1)
i = islice_selected_rec(1)
i = ispec_selected_rec(1)
-
+
endif ! adjoint
! master prints out source time function to file
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_po.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_po.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_po.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -224,11 +224,11 @@
if (SIMULATION_TYPE == 3) then
! to do
stop 'compute_coupling_elastic_po() : adjoint run not implemented yet'
-
+
! dummy to avoid compiler warnings
- iglob = NGLOB_ADJOINT
- iglob = NSPEC_ADJOINT
-
+ iglob = NGLOB_ADJOINT
+ iglob = NSPEC_ADJOINT
+
endif ! adjoint
!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_poroelastic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_poroelastic_el.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_poroelastic_el.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -360,12 +360,12 @@
if (SIMULATION_TYPE == 3) then
! to do
stop 'compute_coupling_poroelastic_el() : adjoint run not implemented yet'
-
+
! dummy to avoid compiler warnings
- iglob = NGLOB_ADJOINT
- iglob = NSPEC_ADJOINT
+ iglob = NGLOB_ADJOINT
+ iglob = NSPEC_ADJOINT
endif ! adjoint
-
+
!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
!!! can merge these loops because NGLLX = NGLLY = NGLLZ do
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -84,10 +84,17 @@
call acoustic_enforce_free_surf_cuda(Mesh_pointer,SIMULATION_TYPE,ABSORB_FREE_SURFACE)
endif
- if(PML) then
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
! enforces free surface on PML elements
- if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ ! note:
+ ! PML routines are not implemented as CUDA kernels, we just transfer the fields
+ ! from the GPU to the CPU and vice versa
+
+ ! transfers potentials to the CPU
+ if(GPU_MODE) call transfer_fields_ac_from_device(NGLOB_AB,potential_acoustic, &
+ potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
+
call PML_acoustic_enforce_free_srfc(NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
ibool,free_surface_ijk,free_surface_ispec, &
@@ -97,11 +104,10 @@
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi2_t_dot_dot,&
chi3_dot_dot,chi4_dot_dot)
- endif
! transfers potentials back to GPU
if(GPU_MODE) call transfer_fields_ac_to_device(NGLOB_AB,potential_acoustic, &
- potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
+ potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
endif
! distinguishes two runs: for elements on MPI interfaces, and elements within the partitions
@@ -191,9 +197,13 @@
num_PML_ispec,PML_ispec,ispec_is_PML_inum,&
chi1_dot,chi2_t,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi3_dot_dot,chi4_dot_dot)
+
+ ! transfers potentials back to GPU
+ if(GPU_MODE) call transfer_fields_ac_to_device(NGLOB_AB,potential_acoustic, &
+ potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
endif
else
- ! Stacey boundary conditions
+ ! Stacey boundary conditions
call compute_stacey_acoustic(NSPEC_AB,NGLOB_AB, &
potential_dot_dot_acoustic,potential_dot_acoustic, &
ibool,ispec_is_inner,phase_is_inner, &
@@ -209,42 +219,42 @@
! elastic coupling
if(ELASTIC_SIMULATION ) then
if( num_coupling_ac_el_faces > 0 ) then
- if( SIMULATION_TYPE == 1 ) then
- ! forward definition: \bfs=\frac{1}{\rho}\bfnabla\phi
- call compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, &
- ibool,displ,potential_dot_dot_acoustic, &
+ if( .NOT. GPU_MODE ) then
+ if( SIMULATION_TYPE == 1 ) then
+ ! forward definition: \bfs=\frac{1}{\rho}\bfnabla\phi
+ call compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, &
+ ibool,displ,potential_dot_dot_acoustic, &
+ num_coupling_ac_el_faces, &
+ coupling_ac_el_ispec,coupling_ac_el_ijk, &
+ coupling_ac_el_normal, &
+ coupling_ac_el_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ else
+ ! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
+ ! adjoint definition: \partial_t^2 \bfs^\dagger=-\frac{1}{\rho}\bfnabla\phi^\dagger
+ call compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, &
+ ibool,-accel_adj_coupling,potential_dot_dot_acoustic, &
+ num_coupling_ac_el_faces, &
+ coupling_ac_el_ispec,coupling_ac_el_ijk, &
+ coupling_ac_el_normal, &
+ coupling_ac_el_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ endif
+ ! adjoint/kernel simulations
+ if( SIMULATION_TYPE == 3 ) &
+ call compute_coupling_acoustic_el(NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ ibool,b_displ,b_potential_dot_dot_acoustic, &
num_coupling_ac_el_faces, &
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
+
else
- ! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
- ! adjoint definition: \partial_t^2 \bfs^\dagger=-\frac{1}{\rho}\bfnabla\phi^\dagger
- call compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, &
- ibool,-accel_adj_coupling,potential_dot_dot_acoustic, &
- num_coupling_ac_el_faces, &
- coupling_ac_el_ispec,coupling_ac_el_ijk, &
- coupling_ac_el_normal, &
- coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner)
- endif
- ! adjoint/kernel simulations
- if( SIMULATION_TYPE == 3 ) &
- call compute_coupling_acoustic_el(NSPEC_ADJOINT,NGLOB_ADJOINT, &
- ibool,b_displ,b_potential_dot_dot_acoustic, &
- num_coupling_ac_el_faces, &
- coupling_ac_el_ispec,coupling_ac_el_ijk, &
- coupling_ac_el_normal, &
- coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner)
- endif
- else
- ! on GPU
- if( num_coupling_ac_el_faces > 0 ) &
+ ! on GPU
call compute_coupling_ac_el_cuda(Mesh_pointer,phase_is_inner, &
num_coupling_ac_el_faces,SIMULATION_TYPE)
-
+ endif ! GPU_MODE
endif
endif
@@ -252,7 +262,7 @@
if(POROELASTIC_SIMULATION ) then
if( num_coupling_ac_po_faces > 0 ) then
if( SIMULATION_TYPE == 1 ) then
- call compute_coupling_acoustic_po(NSPEC_AB,NGLOB_AB, &
+ call compute_coupling_acoustic_po(NSPEC_AB,NGLOB_AB, &
ibool,displs_poroelastic,displw_poroelastic, &
potential_dot_dot_acoustic, &
num_coupling_ac_po_faces, &
@@ -261,10 +271,10 @@
coupling_ac_po_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
else
- stop 'not implemented yet'
+ stop 'not implemented yet'
endif
if( SIMULATION_TYPE == 3 ) &
- stop 'not implemented yet'
+ stop 'not implemented yet'
endif
endif
@@ -443,6 +453,10 @@
! updates potential_dot_acoustic and potential_dot_dot_acoustic inside PML region for plotting seismograms/movies
if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ ! transfers potentials to CPU
+ if(GPU_MODE) call transfer_fields_ac_from_device(NGLOB_AB,potential_acoustic, &
+ potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
+
call PML_acoustic_update_potentials(NGLOB_AB,NSPEC_AB, &
ibool,ispec_is_acoustic, &
potential_dot_acoustic,potential_dot_dot_acoustic,&
@@ -454,7 +468,6 @@
chi1,chi2,chi2_t,chi3,&
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi3_dot_dot,chi4_dot_dot)
- endif
! transfers potentials to GPU
if(GPU_MODE) call transfer_fields_ac_to_device(NGLOB_AB,potential_acoustic, &
@@ -470,15 +483,15 @@
ibool,free_surface_ijk,free_surface_ispec, &
num_free_surface_faces,ispec_is_acoustic)
- if( SIMULATION_TYPE /= 1 ) then
- potential_acoustic_adj_coupling(:) = potential_acoustic(:) &
+ if( SIMULATION_TYPE /= 1 ) then
+ potential_acoustic_adj_coupling(:) = potential_acoustic(:) &
+ deltat * potential_dot_acoustic(:) &
+ deltatsqover2 * potential_dot_dot_acoustic(:)
- endif
+ endif
- ! adjoint simulations
- if (SIMULATION_TYPE == 3) &
- call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT, &
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) &
+ call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT, &
b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
ibool,free_surface_ijk,free_surface_ispec, &
num_free_surface_faces,ispec_is_acoustic)
@@ -489,6 +502,10 @@
if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ ! enforces free surface on PML elements
+ if( GPU_MODE ) call transfer_fields_ac_from_device(NGLOB_AB,potential_acoustic, &
+ potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
+
call PML_acoustic_enforce_free_srfc(NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
ibool,free_surface_ijk,free_surface_ispec, &
@@ -499,7 +516,6 @@
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi2_t_dot_dot,&
chi3_dot_dot,chi4_dot_dot)
- endif
if( GPU_MODE ) call transfer_fields_ac_to_device(NGLOB_AB,potential_acoustic, &
potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -144,14 +144,28 @@
! acoustic coupling
if( ACOUSTIC_SIMULATION ) then
- if( .NOT. GPU_MODE ) then
- call compute_coupling_elastic_ac(NSPEC_AB,NGLOB_AB, &
+ if( num_coupling_ac_el_faces > 0 ) then
+ if( .NOT. GPU_MODE ) then
+ if( SIMULATION_TYPE == 1 ) then
+ ! forward definition: pressure=-potential_dot_dot
+ call compute_coupling_elastic_ac(NSPEC_AB,NGLOB_AB, &
ibool,accel,potential_dot_dot_acoustic, &
num_coupling_ac_el_faces, &
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
+ else
+ ! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
+ ! adoint definition: pressure^\dagger=potential^\dagger
+ call compute_coupling_elastic_ac(NSPEC_AB,NGLOB_AB, &
+ ibool,accel,-potential_acoustic_adj_coupling, &
+ num_coupling_ac_el_faces, &
+ coupling_ac_el_ispec,coupling_ac_el_ijk, &
+ coupling_ac_el_normal, &
+ coupling_ac_el_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ endif
! adjoint simulations
if( SIMULATION_TYPE == 3 ) &
@@ -162,20 +176,38 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
- else
- ! on GPU
- if( num_coupling_ac_el_faces > 0 ) &
+
+ else
+ ! on GPU
call compute_coupling_el_ac_cuda(Mesh_pointer,phase_is_inner, &
- num_coupling_ac_el_faces,SIMULATION_TYPE)
-
- endif
+ num_coupling_ac_el_faces,SIMULATION_TYPE)
+ endif ! GPU_MODE
+ endif ! num_coupling_ac_el_faces
endif
! poroelastic coupling
-! not implemented yet
-! if( POROELASTIC_SIMULATION ) &
-! call compute_coupling_elastic_poro()
+ if( POROELASTIC_SIMULATION ) then
+ call compute_coupling_elastic_po(NSPEC_AB,NGLOB_AB,ibool,&
+ displs_poroelastic,displw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ kappaarraystore,rhoarraystore,mustore, &
+ phistore,tortstore,jacobian,&
+ displ,accel,kappastore, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store, &
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_coupling_el_po_faces, &
+ coupling_el_po_ispec,coupling_po_el_ispec, &
+ coupling_el_po_ijk,coupling_po_el_ijk, &
+ coupling_el_po_normal, &
+ coupling_el_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ endif
! adds source term (single-force/moment-tensor solution)
call compute_add_sources_elastic( NSPEC_AB,NGLOB_AB,accel, &
@@ -852,4 +884,3 @@
end subroutine compute_forces_elastic_Dev_sim3
-
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_fluid.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_fluid.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_fluid.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -229,9 +229,9 @@
stop 'compute_forces_fluid() : adjoint run not implemented yet'
! to avoid compiler warning
l = NGLOB_ADJOINT
- l = NSPEC_ADJOINT
+ l = NSPEC_ADJOINT
endif
-
+
! first double loop over GLL points to compute and store gradients
do l = 1,NGLLX
hp1 = hprime_xx(i,l)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_poroelastic.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_poroelastic.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -399,10 +399,10 @@
if(SIMULATION_TYPE == 3) then
! to do
stop 'compute_continuity_disp_po_el() : adjoint run not implemented yet'
-
+
! dummy to avoid compiler warnings
- i = NGLOB_ADJOINT
- j = NSPEC_ADJOINT
+ i = NGLOB_ADJOINT
+ j = NSPEC_ADJOINT
endif
endif !if(icount(iglob) ==1)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -77,11 +77,14 @@
integer :: ispec,iglob,i,j,k,iface,igll
!integer:: reclen1,reclen2
+ ! checks if anything to do
+ if( num_abs_boundary_faces == 0 ) return
+
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
! reads in absorbing boundary array when first phase is running
if( phase_is_inner .eqv. .false. ) then
- ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
! uses fortran routine
!read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
!if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/create_color_image.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/create_color_image.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/create_color_image.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -33,32 +33,32 @@
! USER PARAMETER
! image data output:
- ! type = 1 : velocity V_x component
- ! type = 2 : velocity V_y component
- ! type = 3 : velocity V_z component
- ! type = 4 : velocity V norm
- integer,parameter:: IMAGE_TYPE = 3
+ ! type = 1 : displ/velocity x-component
+ ! type = 2 : displ/velocity y-component
+ ! type = 3 : displ/velocity z-component
+ ! type = 4 : displ/velocity norm
+ integer,parameter:: IMAGE_TYPE = 3 ! 4
! cross-section surface
! cross-section origin point
- real(kind=CUSTOM_REAL),parameter:: section_xorg = 67000.0
- real(kind=CUSTOM_REAL),parameter:: section_yorg = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_zorg = -1000.0
+ real(kind=CUSTOM_REAL),parameter:: section_xorg = 0.0 ! 67000.0
+ real(kind=CUSTOM_REAL),parameter:: section_yorg = 0.0 ! 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_zorg = -100.0 ! 0.0
! cross-section surface normal
- real(kind=CUSTOM_REAL),parameter:: section_nx = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_ny = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_nz = 1.0
+ real(kind=CUSTOM_REAL),parameter:: section_nx = 0.0 !1.0
+ real(kind=CUSTOM_REAL),parameter:: section_ny = 0.0 !0.0
+ real(kind=CUSTOM_REAL),parameter:: section_nz = 1.0 !0.0
! cross-section (in-plane) horizontal-direction
- real(kind=CUSTOM_REAL),parameter:: section_hdirx = 1.0
- real(kind=CUSTOM_REAL),parameter:: section_hdiry = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_hdirz = 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_hdirx = 1.0 ! 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_hdiry = 0.0 !1.0
+ real(kind=CUSTOM_REAL),parameter:: section_hdirz = 0.0 ! 0.0
! cross-section (in-plane) vertical-direction
- real(kind=CUSTOM_REAL),parameter:: section_vdirx = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_vdiry = 1.0
- real(kind=CUSTOM_REAL),parameter:: section_vdirz = 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_vdirx = 0.0 ! 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_vdiry = 1.0 ! 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_vdirz = 0.0 ! 1.0
! non linear display to enhance small amplitudes in color images
real(kind=CUSTOM_REAL), parameter :: POWER_DISPLAY_COLOR = 0.30_CUSTOM_REAL
@@ -884,7 +884,14 @@
val_vector(:) = displ(:,iglob)
endif
else
- veloc_val(:) = veloc(:,iglob)
+ if( SIMULATION_TYPE == 3 ) then
+ ! to display re-constructed wavefield
+ !val_vector(:) = b_veloc(:,iglob)
+ ! to display adjoint wavefield
+ val_vector(:) = veloc(:,iglob)
+ else
+ val_vector(:) = veloc(:,iglob)
+ endif
endif
! returns with this result
@@ -901,7 +908,15 @@
b_potential_acoustic, val_element,&
hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- ibool,rhostore)
+ ibool,rhostore,GRAVITY)
+ else
+ ! displacement vector
+ call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ potential_acoustic, val_element,&
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ ibool,rhostore,GRAVITY)
+ endif
else
if( SIMULATION_TYPE == 3 ) then
! velocity vector for backward/reconstructed wavefield
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -115,8 +115,8 @@
case( IMODEL_USER_EXTERNAL )
write(IMAIN,'(a)',advance='yes') ' external'
end select
-
- write(IMAIN,*)
+
+ write(IMAIN,*)
endif
! reads in numbers of spectral elements and points for this process' domain
@@ -220,7 +220,7 @@
write(IMAIN,*) 'error: number of MPI processors actually run on: ',sizeprocs
print*
print*, 'error specfem3D: number of processors supposed to run on: ',NPROC
- print*, 'error specfem3D: number of MPI processors actually run on: ',sizeprocs
+ print*, 'error specfem3D: number of MPI processors actually run on: ',sizeprocs
print*
endif
call exit_MPI(myrank,'wrong number of MPI processes')
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -178,20 +178,38 @@
ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
ihours_total,iminutes_total,iseconds_total,int_t_total
+ ! maximum of the norm of the displacement
+ real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all ! elastic
+ real(kind=CUSTOM_REAL) Usolidnormp,Usolidnormp_all ! acoustic
+ real(kind=CUSTOM_REAL) Usolidnorms,Usolidnorms_all ! solid poroelastic
+ real(kind=CUSTOM_REAL) Usolidnormw,Usolidnormw_all ! fluid (w.r.t.s) poroelastic
+
+ ! norm of the backward displacement
+ real(kind=CUSTOM_REAL) b_Usolidnorm, b_Usolidnorm_all
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!chris: Rewrite to get norm for each material when coupled simulations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! compute maximum of norm of displacement in each slice
- if( ELASTIC_SIMULATION ) &
+ if( ELASTIC_SIMULATION ) then
if( GPU_MODE) then
! way 2: just get maximum of field from GPU
call get_norm_elastic_from_device(Usolidnorm,Mesh_pointer,1)
else
Usolidnorm = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2 + displ(3,:)**2))
endif
- if( ACOUSTIC_SIMULATION ) &
+ endif
+
+ if( ACOUSTIC_SIMULATION ) then
+ if(GPU_MODE) then
+ ! way 2: just get maximum of field from GPU
+ call get_norm_acoustic_from_device(Usolidnormp,Mesh_pointer,1)
+ else
Usolidnormp = maxval(abs(potential_dot_dot_acoustic(:)))
+ endif
+ endif
+
if( POROELASTIC_SIMULATION ) then
Usolidnorms = maxval(sqrt(displs_poroelastic(1,:)**2 + displs_poroelastic(2,:)**2 + &
displs_poroelastic(3,:)**2))
@@ -204,13 +222,13 @@
if(Usolidnorm > STABILITY_THRESHOLD .or. Usolidnorm < 0) &
call exit_MPI(myrank,'forward simulation became unstable and blew up')
-! compute the maximum of the maxima for all the slices using an MPI reduction
+ ! compute the maximum of the maxima for all the slices using an MPI reduction
call max_all_cr(Usolidnorm,Usolidnorm_all)
call max_all_cr(Usolidnormp,Usolidnormp_all)
call max_all_cr(Usolidnorms,Usolidnorms_all)
call max_all_cr(Usolidnormw,Usolidnormw_all)
-! adjoint simulations
+ ! adjoint simulations
if( SIMULATION_TYPE == 3 ) then
if( ELASTIC_SIMULATION ) then
! way 2
@@ -239,13 +257,13 @@
call max_all_cr(b_Usolidnorm,b_Usolidnorm_all)
endif
-! user output
+ ! user output
if(myrank == 0) then
write(IMAIN,*) 'Time step # ',it
write(IMAIN,*) 'Time: ',sngl((it-1)*DT-t0),' seconds'
-! elapsed time since beginning of the simulation
+ ! elapsed time since beginning of the simulation
tCPU = wtime() - time_start
int_tCPU = int(tCPU)
ihours = int_tCPU / 3600
@@ -254,10 +272,12 @@
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 = ',sngl(tCPU/dble(it))
- if( ELASTIC_SIMULATION ) &
+ if( ELASTIC_SIMULATION ) then
write(IMAIN,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
- if( ACOUSTIC_SIMULATION ) &
+ endif
+ if( ACOUSTIC_SIMULATION ) then
write(IMAIN,*) 'Max norm pressure P in all slices (Pa) = ',Usolidnormp_all
+ endif
if( POROELASTIC_SIMULATION ) then
write(IMAIN,*) 'Max norm displacement vector Us in all slices (m) = ',Usolidnorms_all
write(IMAIN,*) 'Max norm displacement vector W in all slices (m) = ',Usolidnormw_all
@@ -266,7 +286,7 @@
if (SIMULATION_TYPE == 3) write(IMAIN,*) &
'Max norm U (backward) in all slices = ',b_Usolidnorm_all
-! compute estimated remaining simulation time
+ ! compute estimated remaining simulation time
t_remain = (NSTEP - it) * (tCPU/dble(it))
int_t_remain = int(t_remain)
ihours_remain = int_t_remain / 3600
@@ -278,7 +298,7 @@
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
+ ! compute estimated total simulation time
t_total = t_remain + tCPU
int_t_total = int(t_total)
ihours_total = int_t_total / 3600
@@ -297,7 +317,7 @@
endif
write(IMAIN,*)
-! write time stamp file to give information about progression of simulation
+ ! write time stamp file to give information about progression of simulation
write(outputname,"('/timestamp',i6.6)") it
open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(IOUT,*) 'Time step # ',it
@@ -305,10 +325,12 @@
write(IOUT,*) 'Elapsed time in seconds = ',tCPU
write(IOUT,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
- if( ELASTIC_SIMULATION ) &
+ if( ELASTIC_SIMULATION ) then
write(IOUT,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
- if( ACOUSTIC_SIMULATION ) &
+ endif
+ if( ACOUSTIC_SIMULATION ) then
write(IOUT,*) 'Max norm pressure P in all slices (Pa) = ',Usolidnormp_all
+ endif
if( POROELASTIC_SIMULATION ) then
write(IOUT,*) 'Max norm displacement vector Us in all slices (m) = ',Usolidnorms_all
write(IOUT,*) 'Max norm displacement vector W in all slices (m) = ',Usolidnormw_all
@@ -318,20 +340,19 @@
'Max norm U (backward) in all slices = ',b_Usolidnorm_all
! estimation
write(IOUT,*) 'Time steps done = ',it,' out of ',NSTEP
- write(IOUT,*) 'Time steps remaining = ',NSTEP - it
+ write(IOUT,*) 'Time steps remaining = ',NSTEP - it
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
+ ihours_remain,iminutes_remain,iseconds_remain
write(IOUT,*) 'Estimated total run time in seconds = ',t_total
write(IOUT,"(' Estimated total run time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
ihours_total,iminutes_total,iseconds_total
- write(IOUT,*) 'We have done ',sngl(100.d0*dble(it)/dble(NSTEP)),'% of that'
+ write(IOUT,*) 'We have done ',sngl(100.d0*dble(it)/dble(NSTEP)),'% of that'
close(IOUT)
-
-! check stability of the code, exit if unstable
-! negative values can occur with some compilers when the unstable value is greater
-! than the greatest possible floating-point number of the machine
+ ! check stability of the code, exit if unstable
+ ! negative values can occur with some compilers when the unstable value is greater
+ ! than the greatest possible floating-point number of the machine
if(Usolidnorm_all > STABILITY_THRESHOLD .or. Usolidnorm_all < 0 &
.or. Usolidnormp_all > STABILITY_THRESHOLD .or. Usolidnormp_all < 0 &
.or. Usolidnorms_all > STABILITY_THRESHOLD .or. Usolidnorms_all < 0 &
@@ -408,6 +429,9 @@
! time marching potentials
if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ if( GPU_MODE ) call transfer_fields_ac_from_device(NGLOB_AB,potential_acoustic, &
+ potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
+
call PML_acoustic_time_march(NSPEC_AB,NGLOB_AB,ibool,&
potential_acoustic,potential_dot_acoustic,&
deltat,deltatsqover2,deltatover2,&
@@ -420,8 +444,6 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
my_neighbours_ext_mesh,NPROC,&
ispec_is_acoustic)
- endif
- endif
if( GPU_MODE ) call transfer_fields_ac_to_device(NGLOB_AB,potential_acoustic, &
potential_dot_acoustic, potential_dot_dot_acoustic, Mesh_pointer)
@@ -431,10 +453,19 @@
! updates elastic displacement and velocity
if( ELASTIC_SIMULATION ) then
- displ(:,:) = displ(:,:) + deltat*veloc(:,:) + deltatsqover2*accel(:,:)
- veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
- if( SIMULATION_TYPE /= 1 ) accel_adj_coupling(:,:) = accel(:,:)
- accel(:,:) = 0._CUSTOM_REAL
+
+ if(.NOT. GPU_MODE) then
+ ! on CPU
+ displ(:,:) = displ(:,:) + deltat*veloc(:,:) + deltatsqover2*accel(:,:)
+ veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
+ if( SIMULATION_TYPE /= 1 ) accel_adj_coupling(:,:) = accel(:,:)
+ accel(:,:) = 0._CUSTOM_REAL
+ else
+ ! on GPU
+ ! Includes SIM_TYPE 1 & 3 (for noise tomography)
+ call it_update_displacement_cuda(Mesh_pointer, size(displ), deltat, deltatsqover2,&
+ deltatover2, SIMULATION_TYPE, b_deltat, b_deltatsqover2, b_deltatover2)
+ endif
endif
! updates poroelastic displacements and velocities
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_receivers.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_receivers.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -86,17 +86,17 @@
integer ios
double precision,dimension(1) :: altitude_rec,distmin_ele
- !double precision,dimension(4) :: elevation_node,dist_node
+ !double precision,dimension(4) :: elevation_node,dist_node
double precision,dimension(NPROC) :: distmin_ele_all,elevation_all
real(kind=CUSTOM_REAL) :: xloc,yloc,loc_ele,loc_distmin
-
+
double precision, allocatable, dimension(:) :: x_target,y_target,z_target
-
+
double precision, allocatable, dimension(:) :: horiz_dist
double precision, allocatable, dimension(:) :: x_found,y_found,z_found
double precision dist
-
+
double precision xi,eta,gamma,dx,dy,dz,dxi,deta,dgamma
double precision x,y,z
double precision xix,xiy,xiz
@@ -220,11 +220,11 @@
status='unknown',action='write',iostat=ios)
if( ios /= 0 ) &
call exit_mpi(myrank,'error opening file '//trim(OUTPUT_FILES)//'/output_list_stations.txt')
-
+
do irec=1,nrec
write(IOUT_SU,*) x_found(irec),y_found(irec),z_found(irec)
enddo
-
+
close(IOUT_SU)
deallocate(x_found,y_found,z_found)
endif
@@ -286,7 +286,7 @@
read(1,*,iostat=ios) station_name(irec),network_name(irec), &
stlat(irec),stlon(irec),stele(irec),stbur(irec)
-
+
if (ios /= 0) call exit_mpi(myrank, 'Error reading station file '//trim(rec_filename))
! convert station location to UTM
@@ -315,7 +315,7 @@
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
altitude_rec(1) = loc_ele
distmin_ele(1) = loc_distmin
-
+
! ! set distance to huge initial value
! distmin = HUGEVAL
! if(num_free_surface_faces > 0) then
@@ -385,11 +385,11 @@
! end do
! end do
! end if
-
+
! MPI communications to determine the best slice
call gather_all_dp(distmin_ele,1,distmin_ele_all,1,NPROC)
call gather_all_dp(altitude_rec,1,elevation_all,1,NPROC)
-
+
if(myrank == 0) then
iproc = minloc(distmin_ele_all)
altitude_rec(1) = elevation_all(iproc(1))
@@ -519,7 +519,7 @@
(y_target(irec)>=ymin_ELE .and. y_target(irec)<=ymax_ELE) .and. &
(z_target(irec)>=zmin_ELE .and. z_target(irec)<=zmax_ELE) ) then
! we find the element (ispec) which "may" contain the receiver (irec)
- ! so we only need to compute distances
+ ! so we only need to compute distances
!(which is expensive because of "dsqrt") within those elements
ispec_selected_rec(irec) = ispec
do k = kmin_temp,kmax_temp
@@ -994,11 +994,11 @@
status='unknown',action='write',iostat=ios)
if( ios /= 0 ) &
call exit_mpi(myrank,'error opening file '//trim(OUTPUT_FILES)//'/output_list_stations.txt')
-
+
do irec=1,nrec
write(IOUT_SU,*) x_found(irec),y_found(irec),z_found(irec)
enddo
-
+
close(IOUT_SU)
! stores station infos for later runs
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_source.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/locate_source.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -68,7 +68,7 @@
integer iprocloop
integer i,j,k,ispec,iglob,isource
- integer imin,imax,jmin,jmax,kmin,kmax
+ integer imin,imax,jmin,jmax,kmin,kmax
! integer igll,jgll,kgll,inode,iface,iglob_selected,
! integer iselected,jselected,iface_selected,iadjust,jadjust
integer iproc(1)
@@ -103,13 +103,13 @@
double precision x_target_source,y_target_source,z_target_source
double precision,dimension(1) :: altitude_source,distmin_ele
-
+
double precision,dimension(NPROC) :: distmin_ele_all,elevation_all
! double precision,dimension(4) :: elevation_node,dist_node
real(kind=CUSTOM_REAL) :: xloc,yloc,loc_ele,loc_distmin
-
+
integer islice_selected_source(NSOURCES)
! timer MPI
@@ -220,9 +220,9 @@
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
altitude_source(1) = loc_ele
distmin_ele(1) = loc_distmin
-
-
+
+
! ! set distance to huge initial value
! distmin = HUGEVAL
! if(num_free_surface_faces > 0) then
@@ -295,11 +295,11 @@
! end do
! end if
! distmin_ele(1)= distmin
-
+
! MPI communications to determine the best slice
call gather_all_dp(distmin_ele,1,distmin_ele_all,1,NPROC)
call gather_all_dp(altitude_source,1,elevation_all,1,NPROC)
-
+
if(myrank == 0) then
iproc = minloc(distmin_ele_all)
altitude_source(1) = elevation_all(iproc(1))
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/noise_tomography.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/noise_tomography.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -134,8 +134,16 @@
! output parameters
real(kind=CUSTOM_REAL) :: normal_x_noise_out,normal_y_noise_out,normal_z_noise_out,mask_noise_out
! local parameters
- real(kind=CUSTOM_REAL) :: ldummy
+ !PB VARIABLES TO DEFINE THE REGION OF NOISE
+ real(kind=CUSTOM_REAL) :: xcoord,ycoord,zcoord !,xcoord_center,ycoord_center
+ real :: lon,lat,colat,lon_cn,lat_cn,dsigma,d,dmax
+ ! coordinates "x/y/zcoord_in" actually contain r theta phi, therefore convert back to x y z
+ ! call rthetaphi_2_xyz(xcoord,ycoord,zcoord, xcoord_in,ycoord_in,zcoord_in)
+ xcoord=xcoord_in
+ ycoord=ycoord_in
+ zcoord=zcoord_in
+
!PB NOT UNIF DISTRIBUTION OF NOISE ON THE SURFACE OF A SPHERE
!PB lon lat colat ARE IN RADIANS SINCE ARE OBTAINED FROM CARTESIAN COORDINATES
!PB lon_cn lat_cn (cn = CENTER OF NOISE REGION) IF NOT, MUST BE CONVERTED IN RADIANS
@@ -198,14 +206,9 @@
!******************************** change your noise characteristics above ****************************************
!*****************************************************************************************************************
- ! dummy assign to avoid compiler warnings
- ldummy = xcoord_in
- ldummy = ycoord_in
- ldummy = zcoord_in
+ end subroutine noise_distribution_dir_non_uni
- end subroutine noise_distribution_direction
-
end module user_noise_distribution
@@ -447,20 +450,33 @@
! size of single record
reclen=CUSTOM_REAL*NDIM*NGLLSQUARE*NSPEC_TOP
- ! total file size
- filesize = reclen
- filesize = filesize*NSTEP
- write(outputname,"('/proc',i6.6,'_surface_movie')") myrank
- if (NOISE_TOMOGRAPHY==1) call open_file_abs_w(2,trim(LOCAL_PATH)//trim(outputname), &
- len_trim(trim(LOCAL_PATH)//trim(outputname)), &
- filesize)
- if (NOISE_TOMOGRAPHY==2) call open_file_abs_r(2,trim(LOCAL_PATH)//trim(outputname), &
- len_trim(trim(LOCAL_PATH)//trim(outputname)), &
- filesize)
- if (NOISE_TOMOGRAPHY==3) call open_file_abs_r(2,trim(LOCAL_PATH)//trim(outputname), &
- len_trim(trim(LOCAL_PATH)//trim(outputname)), &
- filesize)
+ ! only open files if there are surface faces in this paritition
+ if(NSPEC_TOP .gt. 0) then
+
+ ! check integer size limit: size of b_reclen_field must fit onto an 4-byte integer
+ if( NSPEC_TOP > 2147483647 / (CUSTOM_REAL * NGLLSQUARE * NDIM) ) then
+ print *,'reclen of noise surface_movie needed exceeds integer 4-byte limit: ',reclen
+ print *,' ',CUSTOM_REAL, NDIM, NGLLSQUARE, NSPEC_TOP
+ print*,'bit size fortran: ',bit_size(NSPEC_TOP)
+ call exit_MPI(myrank,"error NSPEC_TOP integer limit")
+ endif
+
+ ! total file size
+ filesize = reclen
+ filesize = filesize*NSTEP
+
+ write(outputname,"('/proc',i6.6,'_surface_movie')") myrank
+ if (NOISE_TOMOGRAPHY==1) call open_file_abs_w(2,trim(LOCAL_PATH)//trim(outputname), &
+ len_trim(trim(LOCAL_PATH)//trim(outputname)), &
+ filesize)
+ if (NOISE_TOMOGRAPHY==2) call open_file_abs_r(2,trim(LOCAL_PATH)//trim(outputname), &
+ len_trim(trim(LOCAL_PATH)//trim(outputname)), &
+ filesize)
+ if (NOISE_TOMOGRAPHY==3) call open_file_abs_r(2,trim(LOCAL_PATH)//trim(outputname), &
+ len_trim(trim(LOCAL_PATH)//trim(outputname)), &
+ filesize)
+ endif
endif
end subroutine check_parameters_noise
@@ -641,12 +657,9 @@
integer(kind=8) :: Mesh_pointer
logical :: GPU_MODE
- ! loops over surface points
- ! get coordinates of surface mesh and surface displacement
- do iface = 1, num_free_surface_faces
+ ! writes out wavefield at surface
+ if( num_free_surface_faces > 0 ) then
- ispec = free_surface_ispec(iface)
-
if(.NOT. GPU_MODE) then
! loops over surface points
! get coordinates of surface mesh and surface displacement
@@ -729,9 +742,8 @@
! reads in ensemble noise sources at surface
if( num_free_surface_faces > 0 ) then
- ! loops over surface points
- ! puts noise distrubution and direction onto the surface points
- do iface = 1, num_free_surface_faces
+ ! read surface movie
+ call read_abs(2,noise_surface_movie,CUSTOM_REAL*NDIM*NGLLSQUARE*num_free_surface_faces,it)
if(GPU_MODE) then
call noise_read_add_surface_movie_cu(Mesh_pointer, noise_surface_movie,NOISE_TOMOGRAPHY)
@@ -746,13 +758,10 @@
ispec = free_surface_ispec(iface)
- accel(1,iglob) = accel(1,iglob) + eta * mask_noise(ipoin) * normal_x_noise(ipoin) &
- * free_surface_jacobian2Dw(igll,iface)
- accel(2,iglob) = accel(2,iglob) + eta * mask_noise(ipoin) * normal_y_noise(ipoin) &
- * free_surface_jacobian2Dw(igll,iface)
- accel(3,iglob) = accel(3,iglob) + eta * mask_noise(ipoin) * normal_z_noise(ipoin) &
- * free_surface_jacobian2Dw(igll,iface) ! wgllwgll_xy(i,j) * jacobian2D_top(i,j,iface)
- enddo
+ do igll = 1, NGLLSQUARE
+ i = free_surface_ijk(1,igll,iface)
+ j = free_surface_ijk(2,igll,iface)
+ k = free_surface_ijk(3,igll,iface)
ipoin = ipoin + 1
iglob = ibool(i,j,k,ispec)
@@ -823,15 +832,9 @@
integer(kind=8) :: Mesh_pointer
logical :: GPU_MODE
- ! noise source strength kernel
- ! to keep similar structure to other kernels, the source strength kernel is saved as a volumetric kernel
- ! but only updated at the surface, because the noise is generated there
- ipoin = 0
+ ! updates contribution to noise strength kernel
+ if( num_free_surface_faces > 0 ) then
- ! loops over surface points
- ! puts noise distrubution and direction onto the surface points
- do iface = 1, num_free_surface_faces
-
! read surface movie, needed for Sigma_kl
call read_abs(2,noise_surface_movie,CUSTOM_REAL*NDIM*NGLLSQUARE*num_free_surface_faces,it)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.F90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.F90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -139,6 +139,18 @@
if(FIX_UNDERFLOW_PROBLEM) displ(:,:) = VERYSMALLVAL
endif
+ ! initialize poroelastic arrays to zero/verysmallvall
+ if( POROELASTIC_SIMULATION ) then
+ displs_poroelastic(:,:) = 0._CUSTOM_REAL
+ velocs_poroelastic(:,:) = 0._CUSTOM_REAL
+ accels_poroelastic(:,:) = 0._CUSTOM_REAL
+ displw_poroelastic(:,:) = 0._CUSTOM_REAL
+ velocw_poroelastic(:,:) = 0._CUSTOM_REAL
+ accelw_poroelastic(:,:) = 0._CUSTOM_REAL
+ ! put negligible initial value to avoid very slow underflow trapping
+ if(FIX_UNDERFLOW_PROBLEM) displs_poroelastic(:,:) = VERYSMALLVAL
+ if(FIX_UNDERFLOW_PROBLEM) displw_poroelastic(:,:) = VERYSMALLVAL
+ endif
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
@@ -206,7 +218,13 @@
write(IMAIN,*) ' time step: ',sngl(DT),' s'
write(IMAIN,*) 'number of time steps: ',NSTEP
write(IMAIN,*) 'total simulated time: ',sngl(NSTEP*DT),' seconds'
+ write(IMAIN,*) 'start time:',sngl(-t0),' seconds'
write(IMAIN,*)
+
+ !daniel debug: time estimation
+ ! elastic elements: time per element t_per_element = 1.40789368e-05 s
+ ! total time = nspec * nstep * t_per_element
+
endif
! prepares ADJOINT simulations
@@ -269,9 +287,6 @@
endif ! ELASTIC_SIMULATION
if(POROELASTIC_SIMULATION) then
-
- stop 'poroelastic simulation not implemented yet'
- ! but would be something like this...
call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_solid_poroelastic, &
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -176,7 +176,8 @@
if( ier /= 0 ) stop 'error allocating array one_minus_sum_beta etc.'
! reads mass matrices
- read(27) rmass
+ read(27,iostat=ier) rmass
+ if( ier /= 0 ) stop 'error reading in array rmass'
if( OCEANS ) then
! ocean mass matrix
@@ -227,6 +228,9 @@
! poroelastic
call any_all_l( ANY(ispec_is_poroelastic), POROELASTIC_SIMULATION )
if( POROELASTIC_SIMULATION ) then
+
+ if( GPU_MODE ) call exit_mpi(myrank,'POROELASTICITY not supported by GPU mode yet...')
+
! displacement,velocity,acceleration for the solid (s) & fluid (w) phases
allocate(displs_poroelastic(NDIM,NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array displs_poroelastic'
@@ -425,6 +429,42 @@
if(num_phase_ispec_poroelastic > 0 ) read(27) phase_ispec_inner_poroelastic
endif
+! mesh coloring for GPUs
+ if( USE_MESH_COLORING_GPU ) then
+ ! acoustic domain colors
+ if( ACOUSTIC_SIMULATION ) then
+ read(27) num_colors_outer_acoustic,num_colors_inner_acoustic
+
+ allocate(num_elem_colors_acoustic(num_colors_outer_acoustic + num_colors_inner_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_acoustic array'
+
+ read(27) num_elem_colors_acoustic
+ endif
+ ! elastic domain colors
+ if( ELASTIC_SIMULATION ) then
+ read(27) num_colors_outer_elastic,num_colors_inner_elastic
+
+ allocate(num_elem_colors_elastic(num_colors_outer_elastic + num_colors_inner_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_elastic array'
+
+ read(27) num_elem_colors_elastic
+ endif
+ else
+ ! allocates dummy arrays
+ if( ACOUSTIC_SIMULATION ) then
+ num_colors_outer_acoustic = 0
+ num_colors_inner_acoustic = 0
+ allocate(num_elem_colors_acoustic(num_colors_outer_acoustic + num_colors_inner_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_acoustic array'
+ endif
+ if( ELASTIC_SIMULATION ) then
+ num_colors_outer_elastic = 0
+ num_colors_inner_elastic = 0
+ allocate(num_elem_colors_elastic(num_colors_outer_elastic + num_colors_inner_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating num_elem_colors_elastic array'
+ endif
+ endif
+
close(27)
! outputs total element numbers
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -940,12 +940,12 @@
! creates additional receiver and source files
if( SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
- ! extracts receiver locations
- filename = trim(OUTPUT_FILES)//'/sr.vtk'
- filename_new = trim(OUTPUT_FILES)//'/receiver.vtk'
+ ! extracts receiver locations
+ filename = trim(OUTPUT_FILES)//'/sr.vtk'
+ filename_new = trim(OUTPUT_FILES)//'/receiver.vtk'
- ! vtk file for receivers only
- write(system_command, &
+ ! vtk file for receivers only
+ write(system_command, &
"('awk ',a1,'{if(NR<5) print $0;if(NR==5)print ',a1,'POINTS',i6,' float',a1,';if(NR>5+',i6,')print $0}',a1,' < ',a,' > ',a)")&
"'",'"',nrec,'"',NSOURCES,"'",trim(filename),trim(filename_new)
@@ -963,4 +963,26 @@
"('awk ',a1,'{if(NR<5) print $0;if(NR==5)print ',a1,'POINTS',i6,' float',a1,';')") &
"'",'"',NSOURCES,'"'
+ !daniel
+ !print*,'command 1:',trim(system_command1)
+
+ write(system_command2, &
+ "('if(NR>5 && NR <6+',i6,')print $0}END{print ',a,'}',a1,' < ',a,' > ',a)") &
+ NSOURCES,'" "',"'",trim(filename),trim(filename_new)
+
+ !print*,'command 2:',trim(system_command2)
+
+ system_command = trim(system_command1)//trim(system_command2)
+
+ !print*,'command:',trim(system_command)
+!daniel:
+! gfortran
+! call system(trim(system_command),system_command_status)
+! ifort
+! ret = system(trim(system_command))
+
+ endif
+ endif
+
+
end subroutine setup_sources_receivers_VTKfile
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -46,6 +46,29 @@
integer :: NPROC
integer :: NSPEC_AB, NGLOB_AB
+! mesh parameters
+ integer, dimension(:,:,:,:), allocatable :: ibool
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore,ystore,zstore
+
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
+
+! material properties
+ ! isotropic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: kappastore,mustore
+
+! CUDA mesh pointer<->integer wrapper
+ integer(kind=8) :: Mesh_pointer
+
+! Global GPU toggle. Set in Par_file
+ logical :: GPU_MODE
+
+! use integer array to store topography values
+ integer :: NX_TOPO,NY_TOPO
+ !double precision :: ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
+ !character(len=100) :: topo_file
+ integer, dimension(:,:), allocatable :: itopo_bathy
+
! absorbing boundary arrays (for all boundaries) - keeps all infos, allowing for irregular surfaces
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: abs_boundary_normal
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: abs_boundary_jacobian2Dw
@@ -137,7 +160,7 @@
integer :: SIMULATION_TYPE
integer :: NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE
integer :: IMODEL
-
+
double precision :: DT
logical :: ATTENUATION,USE_OLSEN_ATTENUATION, &
@@ -189,13 +212,6 @@
! MPI partition surfaces
logical, dimension(:), allocatable :: ispec_is_inner
-! maximum of the norm of the displacement
- real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all ! elastic
- real(kind=CUSTOM_REAL) Usolidnormp,Usolidnormp_all ! acoustic
- real(kind=CUSTOM_REAL) Usolidnorms,Usolidnorms_all ! solid poroelastic
- real(kind=CUSTOM_REAL) Usolidnormw,Usolidnormw_all ! fluid (w.r.t.s) poroelastic
- integer:: Usolidnorm_index(1)
-
! maximum speed in velocity model
real(kind=CUSTOM_REAL):: model_speed_max
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/write_movie_output.f90 2012-05-12 23:41:07 UTC (rev 20102)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/write_movie_output.f90 2012-05-13 20:21:06 UTC (rev 20103)
@@ -170,7 +170,7 @@
if (USE_HIGHRES_FOR_MOVIES) then
do ipoin = 1, NGLLX*NGLLY
iglob = faces_surface_ext_mesh(ipoin,ispec2D)
-
+
! saves norm of displacement,velocity and acceleration vector
if( ispec_is_elastic(ispec) ) then
! norm of displacement
@@ -310,7 +310,7 @@
logical :: is_done
is_done = .false.
-
+
! loops over all gll points from this element
do k=1,NGLLZ
do j=1,NGLLY
@@ -407,14 +407,6 @@
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
ibool,rhostore,GRAVITY)
endif
- else
- ! velocity vector
- call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
- potential_dot_acoustic, val_element,&
- hprime_xx,hprime_yy,hprime_zz, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- ibool,rhostore)
- endif
endif
if (USE_HIGHRES_FOR_MOVIES) then
@@ -972,14 +964,14 @@
! velocity vector
is_done = .false.
-
+
! loops over all gll points from this element
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
! checks if global point is found
if( iglob == ibool(i,j,k,ispec) ) then
-
+
! horizontal displacement
store_val_ux_external_mesh(ipoin) = max(store_val_ux_external_mesh(ipoin),&
abs(displ_element(1,i,j,k)),abs(displ_element(2,i,j,k)))
Added: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/update_headers_change_word_f90.pl
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/update_headers_change_word_f90.pl (rev 0)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/update_headers_change_word_f90.pl 2012-05-13 20:21:06 UTC (rev 20103)
@@ -0,0 +1,114 @@
+#!/usr/bin/perl
+
+#
+# Script to change the version number in f90 codes
+#
+# Author : Dimitri Komatitsch, EPS - Harvard University, May 1998
+#
+
+#
+# read all f90 and F90 (and *.h) files in the current directory
+# f90 files are supposed to have the extension "*.f90" or "*.F90" or "*.h"
+#
+
+#
+# known bug : does the changes also in constant strings, and not only
+# in the code (which is dangerous, but really easier to program...)
+#
+
+#
+# usage: ./update_headers_change_word_f90.pl
+# run in directory root SPECFEM3D/
+#
+
+
+ at objects = `ls src/*/*.f90 src/*/*.F90 src/*/*.h.in src/*/*.h src/*/*.c src/*/*.cu`;
+
+$nlines_total = 0;
+$nlines_noblank = 0;
+$nlines_nocomment = 0;
+
+foreach $name (@objects) {
+ chop $name;
+
+
+ # change tabs to white spaces
+ if( 1 == 1 ){
+ system("expand -2 < $name > _____tutu01_____");
+ $f90name = $name;
+ print STDOUT "Changing word in file $name ...\n";
+
+ open(FILEF77,"<_____tutu01_____");
+ open(FILEF90,">$f90name");
+
+ # open the source file
+ while($line = <FILEF77>) {
+ chop $line;
+
+ # suppress trailing white spaces and carriage return
+ $line =~ s/\s*$//;
+
+ # change the version number and copyright information
+ # $line =~ s#\(c\) California Institute of Technology and University of Pau, October 2007#\(c\) California Institute of Technology and University of Pau, November 2007#og;
+ # $line =~ s#rmass_sigma#rmass_time_integral_of_sigma#og;
+
+ # write the modified line to the output file
+ print FILEF90 "$line\n";
+
+ }
+
+ close(FILEF77);
+ close(FILEF90);
+ }
+
+ # line count
+ if( 1 == 0 ){
+ print STDOUT "file $name : \n";
+
+ # counts all lines in file
+ $l = `wc -l $name | awk '{print \$1}'`;
+ chomp $l;
+ print " lines = $l \n";
+
+ # without blank lines
+ $lb = 0;
+ # without comments
+ $lc = 0;
+ system("expand -2 < $name > _____tutu01_____");
+ open(FILEF77,"<_____tutu01_____");
+ # open the source file
+ while($line = <FILEF77>) {
+ chop $line;
+ chomp $line;
+ # remove whitespace at start
+ $line =~ s/^\s+//;
+ # remove whitespace at end
+ $line =~ s/\s+$//;
+ # write the modified line to the output file
+ if( $line ne ""){ $lb = $lb + 1; }
+ if( ($line ne "") && (substr($line,0,1) ne "!") && (substr($line,0,1) ne "/") ){ $lc = $lc + 1; }
+ }
+ close(FILEF77);
+ print " lines (no blank) = $lb \n";
+ print " lines (no comment) = $lc \n";
+
+ # summations
+ $nlines_total = $nlines_total + $l;
+ $nlines_noblank = $nlines_noblank + $lb;
+ $nlines_nocomment = $nlines_nocomment + $lc;
+ print " total = $nlines_total \n";
+ print " total (no blank) = $nlines_noblank \n";
+ print " total (no comment) = $nlines_nocomment \n";
+ }
+}
+
+#line count output
+if( 1 == 0 ){
+ print "\ntotal number of lines: \n";
+ print " lines = $nlines_total \n";
+ print " lines (no blank) = $nlines_noblank \n";
+ print " lines (no comment) = $nlines_nocomment \n\n";
+}
+
+system("rm -f _____tutu01_____");
+
Property changes on: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/update_headers_change_word_f90.pl
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list