[cig-commits] r18842 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases shared specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Aug 22 15:12:55 PDT 2011
Author: danielpeter
Date: 2011-08-22 15:12:55 -0700 (Mon, 22 Aug 2011)
New Revision: 18842
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
Log:
fixes bug with large absorption files (>2GB) for kernel simulations
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -389,8 +389,36 @@
call write_VTK_data_elem_i(nspec,nglob, &
xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
v_tmp_i,filename)
+
+ deallocate(iglob_tmp,v_tmp_i)
endif
+ ! saves free surface points
+ if( num_free_surface_faces > 0 ) then
+ ! saves free surface interface points
+ allocate( iglob_tmp(NGLLSQUARE*num_free_surface_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_tmp'
+ inum = 0
+ iglob_tmp(:) = 0
+ do i=1,num_free_surface_faces
+ do j=1,NGLLSQUARE
+ inum = inum+1
+ iglob_tmp(inum) = ibool(free_surface_ijk(1,j,i), &
+ free_surface_ijk(2,j,i), &
+ free_surface_ijk(3,j,i), &
+ free_surface_ispec(i) )
+ enddo
+ enddo
+ filename = prname(1:len_trim(prname))//'free_surface'
+ call write_VTK_data_points(nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iglob_tmp,NGLLSQUARE*num_free_surface_faces, &
+ filename)
+
+ deallocate(iglob_tmp)
+ endif
+
+
!! saves 1. MPI interface
! if( num_interfaces_ext_mesh >= 1 ) then
! filename = prname(1:len_trim(prname))//'MPI_1_points'
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/write_c_binary.c 2011-08-22 22:12:55 UTC (rev 18842)
@@ -26,6 +26,9 @@
!=====================================================================
*/
+// for large files
+#define _FILE_OFFSET_BITS 64
+
// after Brian's function
#include "config.h"
@@ -75,7 +78,7 @@
Jul 18, 2003
- uses functions fopen/fread/fwrite for binary file I/O
-
+
--------------------------------------- */
#define __USE_GNU
@@ -97,7 +100,7 @@
// absorbing files: instead of passing file descriptor, we use the array index
// index 0 for elastic domain file
-// index 1 for acoutic domain file
+// index 1 for acoustic domain file
// (reserved, but unused yet) index 2 - for NOISE_TOMOGRAPHY (SURFACE_MOVIE)
#define ABS_FILEID 3 // number of absorbing files, note: in C we start indexing arrays from 0
@@ -109,7 +112,7 @@
//void
//FC_FUNC_(open_file_abs_r_fbin,OPEN_FILE_ABS_R_FBIN)(int *fid, char *filename,int *length, int *filesize){
-void open_file_abs_r_fbin(int *fid, char *filename,int *length, int *filesize){
+void open_file_abs_r_fbin(int *fid, char *filename,int *length, long long *filesize){
// opens file for read access
@@ -118,13 +121,14 @@
char * fncopy;
char * blank;
FILE *ft;
-
+ int ret;
+
// checks filesize
if( *filesize == 0 ){
perror("Error file size for reading");
exit(EXIT_FAILURE);
}
-
+
// Trim the file name.
fncopy = strndup(filename, *length);
blank = strchr(fncopy, ' ');
@@ -132,14 +136,40 @@
fncopy[blank - fncopy] = '\0';
}
+/*
+//daniel: debug checks file size
+// see:
+//https://www.securecoding.cert.org/confluence/display/seccode/FIO19-C.+Do+not+use+fseek()+and+ftell()+to+compute+the+size+of+a+file
+ printf("file size: %lld \n",*filesize);
+ int fd;
+ struct stat stbuf;
+ long long size;
+ fd = open(fncopy, O_RDONLY);
+ if(fd == -1) {
+ fprintf(stderr, "Error opening file: %s exiting\n", fncopy);
+ exit(-1);
+ }
+ if( fstat(fd, &stbuf) == 0 ){
+ size = stbuf.st_size;
+ printf("file size found is: %lld (Bytes) \n",size);
+ }
+ close(fd);
+*/
+
// opens file
- ft = fopen( fncopy, "r+" );
+ //ft = fopen( fncopy, "r+" );
+ ft = fopen( fncopy, "rb+" ); // read binary file
if( ft == NULL ) { perror("fopen"); exit(-1); }
+
// sets mode for full buffering
work_buffer[*fid] = (char *)malloc(MAX_B);
- setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
-
+ ret = setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
+ if( ret != 0 ){
+ perror("Error setting working buffer");
+ exit(EXIT_FAILURE);
+ }
+
// stores file index id fid: from 0 to 8
fp_abs[*fid] = ft;
@@ -148,7 +178,7 @@
//void
//FC_FUNC_(open_file_abs_w_fbin,OPEN_FILE_ABS_W_FBIN)(int *fid, char *filename, int *length, int *filesize){
-void open_file_abs_w_fbin(int *fid, char *filename, int *length, int *filesize){
+void open_file_abs_w_fbin(int *fid, char *filename, int *length, long long *filesize){
// opens file for write access
@@ -157,10 +187,11 @@
char * fncopy;
char * blank;
FILE *ft;
-
+ int ret;
+
// checks filesize
if( *filesize == 0 ){
- perror("Error file size for reading");
+ perror("Error file size for writing");
exit(EXIT_FAILURE);
}
@@ -172,13 +203,18 @@
}
// opens file
- ft = fopen( fncopy, "w+" );
+ //ft = fopen( fncopy, "w+" );
+ ft = fopen( fncopy, "wb+" ); // write binary file
if( ft == NULL ) { perror("fopen"); exit(-1); }
// sets mode for full buffering
work_buffer[*fid] = (char *)malloc(MAX_B);
- setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
-
+ ret = setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
+ if( ret != 0 ){
+ perror("Error setting working buffer");
+ exit(EXIT_FAILURE);
+ }
+
// stores file index id fid: from 0 to 8
fp_abs[*fid] = ft;
@@ -200,7 +236,7 @@
//void
//FC_FUNC_(write_abs_fbin,WRITE_ABS_FBIN)(int *fid, void *buffer, int *length, int *index){
-void write_abs_fbin(int *fid, void *buffer, int *length, int *index){
+void write_abs_fbin(int *fid, char *buffer, int *length, int *index){
// writes binary file data in chunks of MAX_B
@@ -216,14 +252,38 @@
buf = buffer;
ret = 0;
- //float dat[2];
- //memcpy(dat,buffer,*length);
- //printf("buffer: %f %f\n",dat[0],dat[1]);
-
+/*
+//daniel: debug
+ float dat[*length/4];
+ memcpy(dat,buffer,*length);
+ printf("buffer length: %d %d\n",*length,*index);
+ printf("buffer size: %d %d \n",sizeof(dat),sizeof(buffer));
+ int i;
+ for(i=0;i< 50;i++){
+ printf("buffer: %d %e \n",i,dat[i]);
+ }
+
+ // positions file pointer (for reverse time access)
+ // make sure to use 64-bit arithmetic to avoid overflow for very large files
+ long long pos,cur;
+
+ pos = ((long long)*length) * (*index -1 );
+ cur = ftell(ft);
+
+ printf("current position: %d %lld %lld \n",*fid,cur,pos);
+ ret = fseek(ft, pos , SEEK_SET);
+ if ( ret != 0 ) {
+ perror("Error fseek");
+ exit(EXIT_FAILURE);
+ }
+ */
+
+
// writes items of maximum MAX_B to the file
while (remlen > 0){
itemlen = MIN(remlen,MAX_B);
+ // note: we want to write out exactly *itemlen* bytes
ret = fwrite(buf,1,itemlen,ft);
if (ret > 0){
donelen = donelen + ret;
@@ -235,30 +295,44 @@
}
}
+//daniel: debug
+// printf("buffer done length: %d %d\n",donelen,*length);
+
+
}
//void
//FC_FUNC_(read_abs_fbin,READ_ABS_FBIN)(int *fid, void *buffer, int *length, int *index){
-void read_abs_fbin(int *fid, void *buffer, int *length, int *index){
+void read_abs_fbin(int *fid, char *buffer, int *length, int *index){
// reads binary file data in chunks of MAX_B
FILE *ft;
- int ret,itemlen,remlen,donelen,pos;
+ int ret,itemlen,remlen,donelen;
+ long long pos;
void *buf;
// file pointer
ft = fp_abs[*fid];
// positions file pointer (for reverse time access)
- pos = (*length) * (*index -1 );
- fseek(ft, pos , SEEK_SET);
+ // make sure to use 64-bit arithmetic to avoid overflow for very large files
+ pos = ((long long)*length) * (*index -1 );
+
+ ret = fseek(ft, pos , SEEK_SET);
+ if ( ret != 0 ) {
+ perror("Error fseek");
+ exit(EXIT_FAILURE);
+ }
donelen = 0;
remlen = *length;
buf = buffer;
ret = 0;
+ // cleans buffer
+ //memset( buf,0,remlen);
+
// reads items of maximum MAX_B to the file
while (remlen > 0){
@@ -280,9 +354,20 @@
}
}
- //float dat[2];
- //memcpy(dat,buffer,*length);
- //printf("return buffer: %f %f\n",dat[0],dat[1]);
+/*
+//daniel: debug
+ printf("position: %lld %d %d \n",pos,*length,*index);
+ printf("buffer done length: %d %d\n",donelen,*length);
+ float dat[*length/4];
+ memcpy(dat,buffer,*length);
+ printf("return buffer length: %d %d\n",*length,*index);
+ printf("return buffer size: %d %d \n",sizeof(dat),sizeof(buffer));
+ int i;
+ for(i=0;i< 50;i++){
+ printf("return buffer: %d %e \n",i,dat[i]);
+ }
+*/
+
}
@@ -322,11 +407,11 @@
// file descriptors
static int map_fd_abs[ABS_FILEID];
// file sizes
-static int filesize_abs[ABS_FILEID];
+static long long filesize_abs[ABS_FILEID];
//void
//FC_FUNC_(open_file_abs_w_map,OPEN_FILE_ABS_W_MAP)(int *fid, char *filename, int *length, int *filesize){
-void open_file_abs_w_map(int *fid, char *filename, int *length, int *filesize){
+void open_file_abs_w_map(int *fid, char *filename, int *length, long long *filesize){
// opens file for write access
@@ -342,6 +427,15 @@
exit(EXIT_FAILURE);
}
+ /*
+ // daniel: debug check filesize below or above 2 GB
+ // filesize gives bytes needed; 4-byte integer limited to +- 2,147,483,648 bytes ~ 2 GB
+ float s = *filesize / 1024. / 1024. / 1024.;
+ if( s > 2.0 ){
+ printf("file size bigger than 2 GB: %lld B or %f GB \n",*filesize,s);
+ }
+ */
+
// Trim the file name.
fncopy = strndup(filename, *length);
blank = strchr(fncopy, ' ');
@@ -366,14 +460,13 @@
free(fncopy);
-
/* Stretch the file size to the size of the (mmapped) array of ints
*/
filesize_abs[*fid] = *filesize;
result = lseek(ft, filesize_abs[*fid] - 1, SEEK_SET);
if (result == -1) {
close(ft);
- perror("Error calling fseek() to 'stretch' the file");
+ perror("Error calling lseek() to 'stretch' the file");
exit(EXIT_FAILURE);
}
@@ -414,7 +507,7 @@
//void
//FC_FUNC_(open_file_abs_r_map,OPEN_FILE_ABS_R_MAP)(int *fid, char *filename,int *length, int *filesize){
-void open_file_abs_r_map(int *fid, char *filename,int *length, int *filesize){
+void open_file_abs_r_map(int *fid, char *filename,int *length, long long *filesize){
// opens file for read access
char * fncopy;
@@ -486,12 +579,13 @@
void write_abs_map(int *fid, char *buffer, int *length , int *index){
char *map;
- int offset;
+ long long offset;
map = map_abs[*fid];
// offset in bytes
- offset = (*index -1 ) * (*length) ;
+ // make sure to use 64-bit arithmetic to avoid overflow for very large files
+ offset = ((long long)*index -1 ) * (*length) ;
// copies buffer to map
memcpy( &map[offset], buffer ,*length );
@@ -503,12 +597,13 @@
void read_abs_map(int *fid, char *buffer, int *length , int *index){
char *map;
- int offset;
+ long long offset;
map = map_abs[*fid];
// offset in bytes
- offset = (*index -1 ) * (*length) ;
+ // make sure to use 64-bit arithmetic to avoid overflow for very large files
+ offset = ((long long)*index -1 ) * (*length) ;
// copies map to buffer
memcpy( buffer, &map[offset], *length );
@@ -534,7 +629,7 @@
*/
void
-FC_FUNC_(open_file_abs_w,OPEN_FILE_ABS_W)(int *fid, char *filename,int *length, int *filesize) {
+FC_FUNC_(open_file_abs_w,OPEN_FILE_ABS_W)(int *fid, char *filename,int *length, long long *filesize) {
#ifdef USE_MAP_FUNCTION
open_file_abs_w_map(fid,filename,length,filesize);
@@ -545,7 +640,7 @@
}
void
-FC_FUNC_(open_file_abs_r,OPEN_FILE_ABS_R)(int *fid, char *filename,int *length, int *filesize) {
+FC_FUNC_(open_file_abs_r,OPEN_FILE_ABS_R)(int *fid, char *filename,int *length, long long *filesize) {
#ifdef USE_MAP_FUNCTION
open_file_abs_r_map(fid,filename,length,filesize);
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -80,13 +80,13 @@
logical :: ibool_read_adj_arrays
integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
- real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearray
! local parameters
double precision :: f0
double precision :: stf
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:,:),allocatable:: adj_sourcearray
real(kind=CUSTOM_REAL) stf_used,stf_used_total_all,time_source
- integer :: isource,iglob,ispec,i,j,k
+ integer :: isource,iglob,ispec,i,j,k,ier
integer :: irec_local,irec
! plotting source time function
@@ -231,6 +231,10 @@
! this must be done carefully, otherwise the adjoint sources may be added twice
if (ibool_read_adj_arrays .and. (.not. phase_is_inner)) then
+ ! allocates temporary source array
+ allocate(adj_sourcearray(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adj_sourcearray'
+
irec_local = 0
do irec = 1, nrec
! compute source arrays
@@ -249,6 +253,7 @@
endif
enddo
+ deallocate(adj_sourcearray)
endif ! if(ibool_read_adj_arrays)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -85,13 +85,13 @@
logical :: ibool_read_adj_arrays
integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY
real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
- real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearray
! local parameters
double precision :: f0
double precision :: stf
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:,:),allocatable:: adj_sourcearray
real(kind=CUSTOM_REAL) stf_used,stf_used_total_all,time_source
- integer :: isource,iglob,i,j,k,ispec
+ integer :: isource,iglob,i,j,k,ispec,ier
integer :: irec_local,irec
! plotting source time function
@@ -215,6 +215,10 @@
! this must be done carefully, otherwise the adjoint sources may be added twice
if (ibool_read_adj_arrays .and. (.not. phase_is_inner)) then
+ ! allocates temporary source array
+ allocate(adj_sourcearray(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adj_sourcearray'
+
irec_local = 0
do irec = 1, nrec
! compute source arrays
@@ -233,7 +237,8 @@
endif
enddo
-
+ deallocate(adj_sourcearray)
+
endif ! if(ibool_read_adj_arrays)
if( it < NSTEP ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -122,7 +122,6 @@
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic )
-
if(PML) then
call compute_forces_acoustic_PML(NSPEC_AB,NGLOB_AB, &
ibool,ispec_is_inner,phase_is_inner, &
@@ -303,7 +302,6 @@
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi3_dot_dot,chi4_dot_dot)
-
! enforces free surface (zeroes potentials at free surface)
call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
@@ -329,7 +327,6 @@
chi1_dot_dot,chi2_t_dot_dot,&
chi3_dot_dot,chi4_dot_dot)
-
end subroutine compute_forces_acoustic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -68,22 +68,22 @@
logical:: SAVE_FORWARD
! local parameters
- real(kind=CUSTOM_REAL) :: rhol,cpl,jacobianw
+ real(kind=CUSTOM_REAL) :: rhol,cpl,jacobianw,absorbl
integer :: ispec,iglob,i,j,k,iface,igll
!integer:: reclen1,reclen2
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
- ! reads in absorbing boundary
- ! the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
-
- ! uses fortran routine
- !read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
- !if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
- ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_potential')
- ! uses c routine for faster reading
- call read_abs(1,b_absorb_potential,b_reclen_potential,NSTEP-it+1)
-
+ ! 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
+ ! uses fortran routine
+ !read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
+ !if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_potential')
+ ! uses c routine for faster reading
+ call read_abs(1,b_absorb_potential,b_reclen_potential,NSTEP-it+1)
+ endif
endif !adjoint
! absorbs absorbing-boundary surface using Sommerfeld condition (vanishing field in the outer-space)
@@ -114,16 +114,19 @@
jacobianw = abs_boundary_jacobian2Dw(igll,iface)
! Sommerfeld condition
+ absorbl = potential_dot_acoustic(iglob) * jacobianw / cpl / rhol
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - potential_dot_acoustic(iglob) * jacobianw / cpl / rhol
+ - absorbl
! adjoint simulations
if (SIMULATION_TYPE == 3) then
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- - b_absorb_potential(igll,iface)
+ - b_absorb_potential(igll,iface)
+
else if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
- b_absorb_potential(igll,iface) = potential_dot_acoustic(iglob) * jacobianw / cpl / rhol
+ b_absorb_potential(igll,iface) = absorbl
+
endif !adjoint
enddo
@@ -134,11 +137,13 @@
! adjoint simulations: stores absorbed wavefield part
if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) then
- ! writes out absorbing boundary value
- ! uses fortran routine
- !write(IOABS_AC,rec=it) b_reclen_potential,b_absorb_potential,b_reclen_potential
- ! uses c routine
- call write_abs(1,b_absorb_potential,b_reclen_potential,it)
+ ! writes out absorbing boundary value only when second phase is running
+ if( phase_is_inner .eqv. .true. ) then
+ ! uses fortran routine
+ !write(IOABS_AC,rec=it) b_reclen_potential,b_absorb_potential,b_reclen_potential
+ ! uses c routine
+ call write_abs(1,b_absorb_potential,b_reclen_potential,it)
+ endif
endif
end subroutine compute_stacey_acoustic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -82,16 +82,16 @@
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
- ! reads in absorbing boundary
- ! the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
-
- ! uses fortran routine
- !read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
- !if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
- ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_field')
- ! uses c routine for faster reading
- call read_abs(0,b_absorb_field,b_reclen_field,NSTEP-it+1)
-
+ ! 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
+ ! uses fortran routine
+ !read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
+ !if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_field')
+ ! uses c routine for faster reading
+ call read_abs(0,b_absorb_field,b_reclen_field,NSTEP-it+1)
+ endif
endif !adjoint
@@ -155,11 +155,13 @@
! adjoint simulations: stores absorbed wavefield part
if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) then
- ! writes out absorbing boundary value
- ! uses fortran routine
- !write(IOABS,rec=it) b_reclen_field,b_absorb_field,b_reclen_field
- ! uses c routine
- call write_abs(0,b_absorb_field,b_reclen_field,it)
+ ! writes out absorbing boundary value only when second phase is running
+ if( phase_is_inner .eqv. .true. ) then
+ ! uses fortran routine
+ !write(IOABS,rec=it) b_reclen_field,b_absorb_field,b_reclen_field
+ ! uses c routine
+ call write_abs(0,b_absorb_field,b_reclen_field,it)
+ endif
endif
end subroutine compute_stacey_elastic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -361,8 +361,10 @@
! creating and filling array num_pixel_loc with the positions of each colored
! pixel owned by the local process (useful for parallel jobs)
- allocate(num_pixel_loc(nb_pixel_loc),stat=ier)
- if( ier /= 0 ) stop 'error allocating array num_pixel_loc'
+ if( nb_pixel_loc > 0 ) then
+ allocate(num_pixel_loc(nb_pixel_loc),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array num_pixel_loc'
+ endif
nb_pixel_loc = 0
do i = 1, NX_IMAGE_color
do j = 1, NZ_IMAGE_color
@@ -372,7 +374,11 @@
endif
enddo
enddo
-
+ ! checks if array is allocated
+ if( nb_pixel_loc > 0 ) then
+ if( .not. allocated(num_pixel_loc) ) call exit_MPI(myrank,'error num_pixel_loc allocation')
+ endif
+
! filling array iglob_image_color, containing info on which process owns which pixels.
allocate(nb_pixel_per_proc(0:NPROC-1),stat=ier)
if( ier /= 0 ) stop 'error allocating array nb_pixel_per_proc'
@@ -387,17 +393,20 @@
if( NPROC > 1 ) then
if (myrank == 0) then
do iproc = 1, NPROC-1
- call recv_i(num_pixel_recv(:,iproc),nb_pixel_per_proc(iproc),iproc,42)
-
- ! stores proc number instead where iglob_image_color wouldn't be defined (=-1)
- do k = 1, nb_pixel_per_proc(iproc)
- j = ceiling(real(num_pixel_recv(k,iproc)) / real(NX_IMAGE_color))
- i = num_pixel_recv(k,iproc) - (j-1)*NX_IMAGE_color
- iglob_image_color(i,j) = iproc
- enddo
+ if( nb_pixel_per_proc(iproc) > 0 ) then
+ call recv_i(num_pixel_recv(:,iproc),nb_pixel_per_proc(iproc),iproc,42)
+ ! stores proc number instead where iglob_image_color wouldn't be defined (=-1)
+ do k = 1, nb_pixel_per_proc(iproc)
+ j = ceiling(real(num_pixel_recv(k,iproc)) / real(NX_IMAGE_color))
+ i = num_pixel_recv(k,iproc) - (j-1)*NX_IMAGE_color
+ iglob_image_color(i,j) = iproc
+ enddo
+ endif
enddo
else
- call send_i(num_pixel_loc(:),nb_pixel_loc,0,42)
+ if( nb_pixel_loc > 0 ) then
+ call send_i(num_pixel_loc(:),nb_pixel_loc,0,42)
+ endif
endif
endif
@@ -414,10 +423,12 @@
if( ier /= 0 ) stop 'error allocating array data_pixel_recv'
data_pixel_recv(:) = 0._CUSTOM_REAL
endif
- allocate(data_pixel_send(nb_pixel_loc),stat=ier)
- if(ier /= 0 ) call exit_mpi(myrank,'error allocating image send data')
- data_pixel_send(:) = 0._CUSTOM_REAL
-
+ if( nb_pixel_loc > 0 ) then
+ allocate(data_pixel_send(nb_pixel_loc),stat=ier)
+ if(ier /= 0 ) call exit_mpi(myrank,'error allocating image send data')
+ data_pixel_send(:) = 0._CUSTOM_REAL
+ endif
+
! handles vp background data
call write_PNM_GIF_vp_background()
@@ -463,17 +474,21 @@
! master collects
if (myrank == 0) then
do iproc = 1, NPROC-1
- call recvv_cr(data_pixel_recv(1),nb_pixel_per_proc(iproc),iproc,43)
- ! fills vp display array
- do k = 1, nb_pixel_per_proc(iproc)
- j = ceiling(real(num_pixel_recv(k,iproc)) / real(NX_IMAGE_color))
- i = num_pixel_recv(k,iproc) - (j-1)*NX_IMAGE_color
- image_color_vp_display(i,j) = data_pixel_recv(k)
- enddo
+ if( nb_pixel_per_proc(iproc) > 0 ) then
+ call recvv_cr(data_pixel_recv(1),nb_pixel_per_proc(iproc),iproc,43)
+ ! fills vp display array
+ do k = 1, nb_pixel_per_proc(iproc)
+ j = ceiling(real(num_pixel_recv(k,iproc)) / real(NX_IMAGE_color))
+ i = num_pixel_recv(k,iproc) - (j-1)*NX_IMAGE_color
+ image_color_vp_display(i,j) = data_pixel_recv(k)
+ enddo
+ endif
enddo
else
- ! slave processes send
- call sendv_cr(data_pixel_send,nb_pixel_loc,0,43)
+ if( nb_pixel_loc > 0 ) then
+ ! slave processes send
+ call sendv_cr(data_pixel_send,nb_pixel_loc,0,43)
+ endif
endif
endif
@@ -528,17 +543,21 @@
if (NPROC > 1) then
if (myrank == 0) then
do iproc = 1, NPROC-1
- call recvv_cr(data_pixel_recv(1),nb_pixel_per_proc(iproc),iproc,43)
- ! distributes on image pixels
- do k = 1, nb_pixel_per_proc(iproc)
- j = ceiling(real(num_pixel_recv(k,iproc)) / real(NX_IMAGE_color))
- i = num_pixel_recv(k,iproc) - (j-1)*NX_IMAGE_color
- image_color_data(i,j) = data_pixel_recv(k)
- enddo
+ if( nb_pixel_per_proc(iproc) > 0 ) then
+ call recvv_cr(data_pixel_recv(1),nb_pixel_per_proc(iproc),iproc,43)
+ ! distributes on image pixels
+ do k = 1, nb_pixel_per_proc(iproc)
+ j = ceiling(real(num_pixel_recv(k,iproc)) / real(NX_IMAGE_color))
+ i = num_pixel_recv(k,iproc) - (j-1)*NX_IMAGE_color
+ image_color_data(i,j) = data_pixel_recv(k)
+ enddo
+ endif
enddo
else
- ! slave processes send
- call sendv_cr(data_pixel_send(1),nb_pixel_loc,0,43)
+ if( nb_pixel_loc > 0 ) then
+ ! slave processes send
+ call sendv_cr(data_pixel_send(1),nb_pixel_loc,0,43)
+ endif
endif
endif
@@ -814,11 +833,11 @@
use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM
use specfem_par_acoustic,only: ACOUSTIC_SIMULATION,potential_dot_acoustic,&
- rhostore,ispec_is_acoustic
- use specfem_par_elastic,only: ELASTIC_SIMULATION,veloc,ispec_is_elastic
+ rhostore,ispec_is_acoustic,b_potential_dot_acoustic
+ use specfem_par_elastic,only: ELASTIC_SIMULATION,veloc,ispec_is_elastic,b_veloc
use specfem_par,only: NSPEC_AB,NGLOB_AB,hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- ibool
+ ibool,SIMULATION_TYPE
implicit none
integer,intent(in) :: iglob,ispec
@@ -831,18 +850,34 @@
! returns first element encountered for iglob index
if( ELASTIC_SIMULATION ) then
if( ispec_is_elastic(ispec) ) then
- veloc_val(:) = veloc(:,iglob)
+ if( SIMULATION_TYPE == 3 ) then
+ veloc_val(:) = b_veloc(:,iglob)
+ else
+ veloc_val(:) = veloc(:,iglob)
+ endif
+
+ ! returns with this result
return
endif
endif
if( ACOUSTIC_SIMULATION ) then
if( ispec_is_acoustic(ispec) ) then
- ! velocity vector
- call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ if( SIMULATION_TYPE == 3 ) then
+ ! velocity vector for backward/reconstructed wavefield
+ call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ b_potential_dot_acoustic, veloc_element,&
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ ibool,rhostore)
+ else
+ ! velocity vector
+ call compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
potential_dot_acoustic, veloc_element,&
hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
ibool,rhostore)
+ endif
+
! returns corresponding iglob velocity entry
do k=1,NGLLZ
do j=1,NGLLY
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -259,9 +259,9 @@
logical :: SAVE_FORWARD
! local parameters
integer :: reclen
+ integer(kind=8) :: filesize
character(len=256) :: outputname
-
if (myrank == 0) then
open(unit=IOUT_NOISE,file=trim(OUTPUT_FILES_PATH)//'NOISE_SIMULATION', &
status='unknown',action='write')
@@ -306,14 +306,23 @@
if (NOISE_TOMOGRAPHY/=0) then
! save/read the surface movie using the same c routine as we do for absorbing boundaries (file ID is 2)
- reclen=CUSTOM_REAL*NDIM*NGLLSQUARE*NSPEC_TOP*NSTEP
+
+ ! 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)),reclen)
+ 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)),reclen)
+ 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)),reclen)
+ len_trim(trim(LOCAL_PATH)//trim(outputname)), &
+ filesize)
endif
end subroutine check_parameters_noise
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -114,7 +114,6 @@
! sets up mass matrices
call prepare_timerun_mass_matrices()
-
! initialize acoustic arrays to zero
if( ACOUSTIC_SIMULATION ) then
potential_acoustic(:) = 0._CUSTOM_REAL
@@ -159,6 +158,9 @@
seismograms_a(:,:,:) = 0._CUSTOM_REAL
endif
+ ! synchronize all the processes
+ call sync_all()
+
! prepares attenuation arrays
call prepare_timerun_attenuation()
@@ -416,6 +418,7 @@
implicit none
! local parameters
integer :: ier
+ integer(kind=8) :: filesize
! seismograms
if (nrec_local > 0 .and. SIMULATION_TYPE == 2 ) then
@@ -533,11 +536,16 @@
allocate(b_absorb_field(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array b_absorb_field'
- b_reclen_field = CUSTOM_REAL * (NDIM * NGLLSQUARE * num_abs_boundary_faces)
+ ! size of single record
+ b_reclen_field = CUSTOM_REAL * NDIM * NGLLSQUARE * num_abs_boundary_faces
+
+ ! total file size
+ filesize = b_reclen_field
+ filesize = filesize*NSTEP
if (SIMULATION_TYPE == 3) then
! opens existing files
-
+
! uses fortran routines for reading
!open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='old',&
! action='read',form='unformatted',access='direct', &
@@ -546,11 +554,10 @@
! uses c routines for faster reading
call open_file_abs_r(0,trim(prname)//'absorb_field.bin', &
len_trim(trim(prname)//'absorb_field.bin'), &
- b_reclen_field*NSTEP)
+ filesize)
else
! opens new file
-
! uses fortran routines for writing
!open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='unknown',&
! form='unformatted',access='direct',&
@@ -559,7 +566,7 @@
! uses c routines for faster writing (file index 0 for acoutic domain file)
call open_file_abs_w(0,trim(prname)//'absorb_field.bin', &
len_trim(trim(prname)//'absorb_field.bin'), &
- b_reclen_field*NSTEP)
+ filesize)
endif
endif
@@ -570,8 +577,21 @@
allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array b_absorb_potential'
- b_reclen_potential = CUSTOM_REAL * (NGLLSQUARE * num_abs_boundary_faces)
+ ! size of single record
+ b_reclen_potential = CUSTOM_REAL * NGLLSQUARE * num_abs_boundary_faces
+
+ ! total file size (two lines to implicitly convert to 8-byte integers)
+ filesize = b_reclen_potential
+ filesize = filesize*NSTEP
+ ! daniel: debug check size limit
+ !if( NSTEP > 2147483648 / b_reclen_potential ) then
+ ! print *,'file size needed exceeds integer 4-byte limit: ',b_reclen_potential,NSTEP
+ ! print *,' ',CUSTOM_REAL, NGLLSQUARE, num_abs_boundary_faces,NSTEP
+ ! print*,'file size fortran: ',filesize
+ ! print*,'file bit size fortran: ',bit_size(filesize)
+ !endif
+
if (SIMULATION_TYPE == 3) then
! opens existing files
! uses fortran routines for reading
@@ -579,10 +599,11 @@
! action='read',form='unformatted',access='direct', &
! recl=b_reclen_potential+2*4,iostat=ier )
!if( ier /= 0 ) call exit_mpi(myrank,'error opening proc***_absorb_potential.bin file')
+
! uses c routines for faster reading
call open_file_abs_r(1,trim(prname)//'absorb_potential.bin', &
len_trim(trim(prname)//'absorb_potential.bin'), &
- b_reclen_potential*NSTEP)
+ filesize)
else
! opens new file
@@ -594,7 +615,7 @@
! uses c routines for faster writing (file index 1 for acoutic domain file)
call open_file_abs_w(1,trim(prname)//'absorb_potential.bin', &
len_trim(trim(prname)//'absorb_potential.bin'), &
- b_reclen_potential*NSTEP)
+ filesize)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-08-20 01:10:07 UTC (rev 18841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-08-22 22:12:55 UTC (rev 18842)
@@ -635,7 +635,9 @@
else
! allocate dummy array in order to be able to use it as a subroutine argument, even if unused
- allocate(adj_sourcearrays(1,1,1,1,1,1),stat=ier)
+ nadj_rec_local = 0
+ NTSTEP_BETWEEN_READ_ADJSRC = 0
+ allocate(adj_sourcearrays(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating dummy array adj_sourcearrays'
endif
More information about the CIG-COMMITS
mailing list