[cig-commits] r22994 - in seismo/3D/SPECFEM3D_GLOBE/trunk: setup src/auxiliaries src/shared src/specfem3D utils/Visualization/VTK_ParaView utils/doubling_brick
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Jan 24 08:03:05 PST 2014
Author: danielpeter
Date: 2014-01-24 08:03:04 -0800 (Fri, 24 Jan 2014)
New Revision: 22994
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/file_io_threads.c
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_adjoint_sources.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Visualization/VTK_ParaView/create_slice_VTK.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v36.h
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v40.h
Log:
uses asynchronuous reading of adjoint sources by use of a separate posix thread, trying to overlap file i/o with kernel computations; adds flag IO_ASYNC_COPY to constants.h.in; adds new files file_io_threads.c and read_adjoint_sources.f90; renames IOVTK to IOUT_VTK
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2014-01-24 16:03:04 UTC (rev 22994)
@@ -58,17 +58,24 @@
logical, parameter :: LOCAL_PATH_IS_ALSO_GLOBAL = . at LOCAL_PATH_IS_ALSO_GLOBAL@.
! input, output and main MPI I/O files
+! note: careful with these unit numbers, we mostly use units in the 40-50 range.
+! cray fortran e.g. reserves 0,5,6 (standard error,input,output units) and 100-102
integer, parameter :: ISTANDARD_OUTPUT = 6
- integer, parameter :: IIN = 40,IOUT = 41,IOUT_SAC = 903
- integer, parameter :: IIN_NOISE = 43,IOUT_NOISE = 44
-! local file unit for output of buffers
- integer, parameter :: IOUT_BUFFERS = 35
+ integer, parameter :: IIN = 40,IOUT = 41
! uncomment this to write messages to a text file
integer, parameter :: IMAIN = 42
! uncomment this to write messages to the screen (slows down the code)
! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
+! I/O unit for noise data files
+ integer, parameter :: IIN_NOISE = 43,IOUT_NOISE = 44
+! I/O unit for adjoint source files
+ integer, parameter :: IIN_ADJ = 45
+! local file unit for output of buffers
+ integer, parameter :: IOUT_BUFFERS = 46
! I/O unit for source and receiver vtk file
- integer, parameter :: IOVTK = 98
+ integer, parameter :: IOUT_VTK = 47
+! I/O unit for sac files
+ integer, parameter :: IOUT_SAC = 48
! I/O Bluegene serial access
logical, parameter :: IO_BLUEGENE_SERIAL = .false.
@@ -281,6 +288,9 @@
!!
!!-----------------------------------------------------------
+! asynchronuous reading of adjoint sources
+ logical, parameter :: IO_ASYNC_COPY = .true.
+
! regular kernel parameters
character(len=*), parameter :: PATHNAME_KL_REG = 'DATA/kl_reg_grid.txt'
integer, parameter :: NM_KL_REG_LAYER = 100
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.F90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.F90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -610,42 +610,42 @@
! VTK
! opens unstructured grid file
write(mesh_file,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'.vtk'
- open(IOVTK,file=mesh_file(1:len_trim(mesh_file)),status='unknown',iostat=ios)
+ open(IOUT_VTK,file=mesh_file(1:len_trim(mesh_file)),status='unknown',iostat=ios)
if( ios /= 0 ) stop 'error opening vtk output file'
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
! points
- write(IOVTK, '(a,i16,a)') 'POINTS ', np, ' float'
+ write(IOUT_VTK, '(a,i16,a)') 'POINTS ', np, ' float'
do i = 1,np
- write(IOVTK,'(3e18.6)') total_dat_xyz(1,i),total_dat_xyz(2,i),total_dat_xyz(3,i)
+ write(IOUT_VTK,'(3e18.6)') total_dat_xyz(1,i),total_dat_xyz(2,i),total_dat_xyz(3,i)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! cells
! note: indices for vtk start at 0
- write(IOVTK,'(a,i12,i12)') "CELLS ",ne,ne*9
+ write(IOUT_VTK,'(a,i12,i12)') "CELLS ",ne,ne*9
do i = 1,ne
- write(IOVTK,'(9i12)') 8,total_dat_con(1,i),total_dat_con(2,i),total_dat_con(3,i),total_dat_con(4,i), &
+ write(IOUT_VTK,'(9i12)') 8,total_dat_con(1,i),total_dat_con(2,i),total_dat_con(3,i),total_dat_con(4,i), &
total_dat_con(5,i),total_dat_con(6,i),total_dat_con(7,i),total_dat_con(8,i)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! VTK
! type: hexahedrons
- write(IOVTK,'(a,i12)') "CELL_TYPES ",ne
- write(IOVTK,'(6i12)') (12,it=1,ne)
- write(IOVTK,*) ""
+ write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",ne
+ write(IOUT_VTK,'(6i12)') (12,it=1,ne)
+ write(IOUT_VTK,*) ""
- write(IOVTK,'(a,i12)') "POINT_DATA ",np
- write(IOVTK,'(a)') "SCALARS "//trim(filename)//" float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "POINT_DATA ",np
+ write(IOUT_VTK,'(a)') "SCALARS "//trim(filename)//" float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do i = 1,np
- write(IOVTK,*) total_dat(i)
+ write(IOUT_VTK,*) total_dat(i)
enddo
- write(IOVTK,*) ""
- close(IOVTK)
+ write(IOUT_VTK,*) ""
+ close(IOUT_VTK)
! free arrays for this region
deallocate(total_dat,total_dat_xyz,total_dat_con)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_VTK_file.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -56,26 +56,26 @@
!write(IMAIN,*) ' vtk file: '
!write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(IOVTK, '(a,i12,a)') 'POINTS ', num_points_globalindices, ' float'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', num_points_globalindices, ' float'
do i=1,num_points_globalindices
iglob = points_globalindices(i)
if( iglob <= 0 .or. iglob > nglob ) then
print*,'error: '//prname_file(1:len_trim(prname_file))//'.vtk'
print*,'error global index: ',iglob,i
- close(IOVTK)
+ close(IOUT_VTK)
stop 'error vtk points file'
endif
- write(IOVTK,'(3e18.6)') xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
- close(IOVTK)
+ close(IOUT_VTK)
end subroutine write_VTK_data_points
@@ -113,27 +113,27 @@
!write(IMAIN,*) ' vtk file: '
!write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
do iglob=1,nglob
- write(IOVTK,*) xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob)
+ write(IOUT_VTK,*) xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! writes out gll-data (velocity) for each element point
- write(IOVTK,'(a,i12)') "POINT_DATA ",nglob
- write(IOVTK,'(a)') "SCALARS glob_data float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "POINT_DATA ",nglob
+ write(IOUT_VTK,'(a)') "SCALARS glob_data float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do iglob=1,nglob
- write(IOVTK,*) glob_values(iglob)
+ write(IOUT_VTK,*) glob_values(iglob)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
- close(IOVTK)
+ close(IOUT_VTK)
end subroutine write_VTK_glob_points
@@ -172,42 +172,42 @@
!write(IMAIN,*) ' vtk file: '
!write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
do i=1,nglob
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! note: indices for vtk start at 0
- write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
do ispec=1,nspec
- write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
+ write(IOUT_VTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
ibool(1,1,NGLLZ,ispec)-1,ibool(NGLLX,1,NGLLZ,ispec)-1,ibool(NGLLX,NGLLY,NGLLZ,ispec)-1,ibool(1,NGLLY,NGLLZ,ispec)-1
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! type: hexahedrons
- write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
- write(IOVTK,'(6i12)') (12,ispec=1,nspec)
- write(IOVTK,*) ""
+ write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOUT_VTK,'(6i12)') (12,ispec=1,nspec)
+ write(IOUT_VTK,*) ""
- write(IOVTK,'(a,i12)') "CELL_DATA ",nspec
- write(IOVTK,'(a)') "SCALARS elem_flag integer"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "CELL_DATA ",nspec
+ write(IOUT_VTK,'(a)') "SCALARS elem_flag integer"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do ispec = 1,nspec
if( elem_flag(ispec) .eqv. .true. ) then
- write(IOVTK,*) 1
+ write(IOUT_VTK,*) 1
else
- write(IOVTK,*) 0
+ write(IOUT_VTK,*) 0
endif
enddo
- write(IOVTK,*) ""
- close(IOVTK)
+ write(IOUT_VTK,*) ""
+ close(IOUT_VTK)
end subroutine write_VTK_data_elem_l
@@ -247,38 +247,38 @@
!write(IMAIN,*) ' vtk file: '
!write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
do i=1,nglob
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! note: indices for vtk start at 0
- write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
do ispec=1,nspec
- write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
+ write(IOUT_VTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
ibool(1,1,NGLLZ,ispec)-1,ibool(NGLLX,1,NGLLZ,ispec)-1,ibool(NGLLX,NGLLY,NGLLZ,ispec)-1,ibool(1,NGLLY,NGLLZ,ispec)-1
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! type: hexahedrons
- write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
- write(IOVTK,'(6i12)') (12,ispec=1,nspec)
- write(IOVTK,*) ""
+ write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOUT_VTK,'(6i12)') (12,ispec=1,nspec)
+ write(IOUT_VTK,*) ""
- write(IOVTK,'(a,i12)') "CELL_DATA ",nspec
- write(IOVTK,'(a)') "SCALARS elem_val integer"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "CELL_DATA ",nspec
+ write(IOUT_VTK,'(a)') "SCALARS elem_val integer"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do ispec = 1,nspec
- write(IOVTK,*) elem_flag(ispec)
+ write(IOUT_VTK,*) elem_flag(ispec)
enddo
- write(IOVTK,*) ""
- close(IOVTK)
+ write(IOUT_VTK,*) ""
+ close(IOUT_VTK)
end subroutine write_VTK_data_elem_i
@@ -317,12 +317,12 @@
real(kind=CUSTOM_REAL) :: rval,thetaval,phival,xval,yval,zval
! write source and receiver VTK files for Paraview
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
do i=1,nglob
!x,y,z store have been converted to r theta phi already, need to revert back for xyz output
@@ -331,21 +331,21 @@
phival = zstore_dummy(i)
call rthetaphi_2_xyz(xval,yval,zval,rval,thetaval,phival)
- !write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
- write(IOVTK,'(3e18.6)') xval,yval,zval
+ !write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xval,yval,zval
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! defines cell on coarse corner points
! note: indices for vtk start at 0
- write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
do ispec=1,nspec
! specific to inner core elements
! exclude fictitious elements in central cube
if(idoubling(ispec) /= IFLAG_IN_FICTITIOUS_CUBE) then
! valid cell
- write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1, &
+ write(IOUT_VTK,'(9i12)') 8,ibool(1,1,1,ispec)-1, &
ibool(NGLLX,1,1,ispec)-1, &
ibool(NGLLX,NGLLY,1,ispec)-1, &
ibool(1,NGLLY,1,ispec)-1, &
@@ -356,7 +356,7 @@
else
! fictitious elements in central cube
! maps cell onto a randomly chosen point
- write(IOVTK,'(9i12)') 8,ibool(1,1,1,1)-1, &
+ write(IOUT_VTK,'(9i12)') 8,ibool(1,1,1,1)-1, &
ibool(1,1,1,1)-1, &
ibool(1,1,1,1)-1, &
ibool(1,1,1,1)-1, &
@@ -367,43 +367,43 @@
endif
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! type: hexahedrons
- write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
- write(IOVTK,'(6i12)') (12,ispec=1,nspec)
- write(IOVTK,*) ""
+ write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOUT_VTK,'(6i12)') (12,ispec=1,nspec)
+ write(IOUT_VTK,*) ""
! x components
- write(IOVTK,'(a,i12)') "POINT_DATA ",nglob
- write(IOVTK,'(a)') "SCALARS x_comp float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "POINT_DATA ",nglob
+ write(IOUT_VTK,'(a)') "SCALARS x_comp float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do i = 1,nglob
- write(IOVTK,*) glob_data(1,i)
+ write(IOUT_VTK,*) glob_data(1,i)
enddo
! y components
- write(IOVTK,'(a)') "SCALARS y_comp float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a)') "SCALARS y_comp float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do i = 1,nglob
- write(IOVTK,*) glob_data(2,i)
+ write(IOUT_VTK,*) glob_data(2,i)
enddo
! z components
- write(IOVTK,'(a)') "SCALARS z_comp float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a)') "SCALARS z_comp float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do i = 1,nglob
- write(IOVTK,*) glob_data(3,i)
+ write(IOUT_VTK,*) glob_data(3,i)
enddo
! norm
- write(IOVTK,'(a)') "SCALARS norm float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a)') "SCALARS norm float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do i = 1,nglob
- write(IOVTK,*) sqrt( glob_data(1,i)*glob_data(1,i) &
+ write(IOUT_VTK,*) sqrt( glob_data(1,i)*glob_data(1,i) &
+ glob_data(2,i)*glob_data(2,i) &
+ glob_data(3,i)*glob_data(3,i))
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
- close(IOVTK)
+ close(IOUT_VTK)
end subroutine write_VTK_data_cr
@@ -499,12 +499,12 @@
if( myrank == 0 ) then
! write source and receiver VTK files for Paraview
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(IOVTK, '(a,i12,a)') 'POINTS ', nglob*NPROCTOT, ' float'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob*NPROCTOT, ' float'
do iproc=0, NPROCTOT-1
do i=1,nglob
@@ -514,15 +514,15 @@
phival = store_val_z_all(i,iproc)
call rthetaphi_2_xyz(xval,yval,zval,rval,thetaval,phival)
- !write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
- write(IOVTK,'(3e18.6)') xval,yval,zval
+ !write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xval,yval,zval
enddo
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! defines cell on coarse corner points
! note: indices for vtk start at 0
- write(IOVTK,'(a,i12,i12)') "CELLS ",nspec*NPROCTOT,nspec*NPROCTOT*9
+ write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec*NPROCTOT,nspec*NPROCTOT*9
do iproc=0, NPROCTOT-1
do ispec=1,nspec
@@ -534,7 +534,7 @@
if(idoubling_all(ispec,iproc) /= IFLAG_IN_FICTITIOUS_CUBE) then
! valid cell
! cell corner ids
- write(IOVTK,'(9i12)') 8,ibool_all(1,1,1,ispec,iproc)-1+iproc*nglob, &
+ write(IOUT_VTK,'(9i12)') 8,ibool_all(1,1,1,ispec,iproc)-1+iproc*nglob, &
ibool_all(NGLLX,1,1,ispec,iproc)-1+iproc*nglob, &
ibool_all(NGLLX,NGLLY,1,ispec,iproc)-1+iproc*nglob, &
ibool_all(1,NGLLY,1,ispec,iproc)-1+iproc*nglob, &
@@ -545,7 +545,7 @@
else
! fictitious elements in central cube
! maps cell onto a randomly chosen point
- write(IOVTK,'(9i12)') 8,ibool_all(1,1,1,1,iproc)-1, &
+ write(IOUT_VTK,'(9i12)') 8,ibool_all(1,1,1,1,iproc)-1, &
ibool_all(1,1,1,1,iproc)-1, &
ibool_all(1,1,1,1,iproc)-1, &
ibool_all(1,1,1,1,iproc)-1, &
@@ -557,51 +557,51 @@
enddo
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! type: hexahedrons
- write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec*NPROCTOT
- write(IOVTK,'(6i12)') (12,ispec=1,nspec*NPROCTOT)
- write(IOVTK,*) ""
+ write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",nspec*NPROCTOT
+ write(IOUT_VTK,'(6i12)') (12,ispec=1,nspec*NPROCTOT)
+ write(IOUT_VTK,*) ""
! x components
- write(IOVTK,'(a,i12)') "POINT_DATA ",nglob*NPROCTOT
- write(IOVTK,'(a)') "SCALARS x_comp float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "POINT_DATA ",nglob*NPROCTOT
+ write(IOUT_VTK,'(a)') "SCALARS x_comp float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do iproc=0, NPROCTOT-1
do i = 1,nglob
- write(IOVTK,*) store_val_ux_all(i,iproc)
+ write(IOUT_VTK,*) store_val_ux_all(i,iproc)
enddo
enddo
! y components
- write(IOVTK,'(a)') "SCALARS y_comp float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a)') "SCALARS y_comp float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do iproc=0, NPROCTOT-1
do i = 1,nglob
- write(IOVTK,*) store_val_uy_all(i,iproc)
+ write(IOUT_VTK,*) store_val_uy_all(i,iproc)
enddo
enddo
! z components
- write(IOVTK,'(a)') "SCALARS z_comp float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a)') "SCALARS z_comp float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do iproc=0, NPROCTOT-1
do i = 1,nglob
- write(IOVTK,*) store_val_uz_all(i,iproc)
+ write(IOUT_VTK,*) store_val_uz_all(i,iproc)
enddo
enddo
! norm
- write(IOVTK,'(a)') "SCALARS norm float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a)') "SCALARS norm float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do iproc=0, NPROCTOT-1
do i = 1,nglob
- write(IOVTK,*) sqrt( store_val_ux_all(i,iproc)**2 &
+ write(IOUT_VTK,*) sqrt( store_val_ux_all(i,iproc)**2 &
+ store_val_uy_all(i,iproc)**2 &
+ store_val_uz_all(i,iproc)**2 )
enddo
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
- close(IOVTK)
+ close(IOUT_VTK)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -151,62 +151,19 @@
implicit none
! local parameters
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: tmp_sourcearray
- integer :: irec,irec_local,i,j,k,iglob,it_sub_adj,itime,ier
+ integer :: irec,irec_local,i,j,k,iglob
integer :: ivec_index
- character(len=150) :: adj_source_file
logical :: ibool_read_adj_arrays
- ! note: we check if nadj_rec_local > 0 before calling this routine
+ ! note: we check if nadj_rec_local > 0 before calling this routine, but better be safe...
+ if( nadj_rec_local == 0 ) return
! figure out if we need to read in a chunk of the adjoint source at this timestep
- it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) ) !chunk_number
- ibool_read_adj_arrays = (((it == it_begin) .or. (mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC) == 0)) &
- .and. (nadj_rec_local > 0))
+ ibool_read_adj_arrays = ( (it == it_begin) .or. (mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC) == 0) )
! needs to read in a new chunk/block of the adjoint source
if(ibool_read_adj_arrays) then
- ! allocates temporary source array
- allocate(tmp_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NTSTEP_BETWEEN_READ_ADJSRC),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating array tmp_sourcearray')
-
- tmp_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
- irec_local = 0
-
- do irec = 1, nrec
- ! check that the source slice number is okay
- if(islice_selected_rec(irec) < 0 .or. islice_selected_rec(irec) > NPROCTOT_VAL-1) then
- if(islice_selected_rec(irec) < 0) call exit_MPI(myrank,'islice < 0')
- if(islice_selected_rec(irec) > NPROCTOT_VAL-1) call exit_MPI(myrank,'islice > NPROCTOT_VAL-1')
- call exit_MPI(myrank,'now: something is wrong with the source slice number in adjoint simulation')
- endif
- ! 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_source_adjoint(myrank,adj_source_file, &
- xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- nu(:,:,irec),tmp_sourcearray, xigll,yigll,zigll,iadjsrc_len(it_sub_adj), &
- iadjsrc,it_sub_adj,NSTEP_SUB_ADJ,NTSTEP_BETWEEN_READ_ADJSRC,DT)
-
- ! stores source array
- ! note: the adj_sourcearrays has a time stepping from 1 to NTSTEP_BETWEEN_READ_ADJSRC
- ! this gets overwritten every time a new block/chunk is read in
- do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
- adj_sourcearrays(:,:,:,:,irec_local,itime) = tmp_sourcearray(:,:,:,:,itime)
- enddo
-
- endif
- enddo
-
- ! checks that number of read sources is valid
- if(irec_local /= nadj_rec_local) &
- call exit_MPI(myrank,'irec_local /= nadj_rec_local in adjoint simulation')
-
- ! frees temporary array
- deallocate(tmp_sourcearray)
+ call read_adjoint_sources()
endif
! adds adjoint sources
@@ -343,7 +300,6 @@
endif
endif
-
endif
end subroutine compute_add_sources_adjoint
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -128,11 +128,11 @@
!================================================================
subroutine compute_arrays_source_adjoint(myrank, adj_source_file, &
- xi_receiver,eta_receiver,gamma_receiver, nu,adj_sourcearray, &
- xigll,yigll,zigll,NSTEP_BLOCK,iadjsrc,it_sub_adj,NSTEP_SUB_ADJ, &
- NTSTEP_BETWEEN_READ_ADJSRC,DT)
+ xi_receiver,eta_receiver,gamma_receiver, nu,adj_sourcearray, &
+ xigll,yigll,zigll,NSTEP_BLOCK,iadjsrc,it_sub_adj,NSTEP_SUB_ADJ, &
+ NTSTEP_BETWEEN_READ_ADJSRC,DT)
- use constants
+ use constants,only: CUSTOM_REAL,SIZE_REAL,NDIM,NGLLX,NGLLY,NGLLZ,IIN_ADJ,R_EARTH
implicit none
@@ -211,7 +211,7 @@
! opens adjoint component file
filename = 'SEM/'//trim(adj_source_file) // '.'// comp(icomp) // '.adj'
- open(unit=IIN,file=trim(filename),status='old',action='read',iostat=ios)
+ open(unit=IIN_ADJ,file=trim(filename),status='old',action='read',iostat=ios)
! note: adjoint source files must be available for all three components E/N/Z, even
! if a component is just zeroed out
@@ -221,11 +221,11 @@
call exit_MPI(myrank,&
'file '//trim(filename)//' not found, please check with your STATIONS_ADJOINT file')
endif
- !if (ios /= 0) cycle ! cycles to next file
+ !if (ios /= 0) cycle ! cycles to next file - this is too error prone and users might easily end up with wrong results
! jumps over unused trace length
do itime =1,it_start-1
- read(IIN,*,iostat=ios) junk,junk
+ read(IIN_ADJ,*,iostat=ios) junk,junk
if( ios /= 0) &
call exit_MPI(myrank,&
'file '//trim(filename)//' has wrong length, please check with your simulation duration')
@@ -244,15 +244,15 @@
endif
! reads in adjoint source trace
- !read(IIN,*,iostat=ios) junk, adj_src(icomp,itime-it_start+1)
- read(IIN,*,iostat=ios) junk, adj_src(icomp,index_i)
+ !read(IIN_ADJ,*,iostat=ios) junk, adj_src(icomp,itime-it_start+1)
+ read(IIN_ADJ,*,iostat=ios) junk, adj_src(icomp,index_i)
if( ios /= 0) &
call exit_MPI(myrank, &
'file '//trim(filename)//' has wrong length, please check with your simulation duration')
enddo
- close(IIN)
+ close(IIN_ADJ)
enddo
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/file_io_threads.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/file_io_threads.c (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/file_io_threads.c 2014-01-24 16:03:04 UTC (rev 22994)
@@ -0,0 +1,420 @@
+/*
+ !=====================================================================
+ !
+ ! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
+ ! --------------------------------------------------
+ !
+ ! Main authors: Dimitri Komatitsch and Jeroen Tromp
+ ! Princeton University, USA
+ ! and CNRS / INRIA / University of Pau, France
+ ! (c) Princeton University and CNRS / INRIA / University of Pau
+ ! August 2013
+ !
+ ! This program is free software; you can redistribute it and/or modify
+ ! it under the terms of the GNU General Public License as published by
+ ! the Free Software Foundation; either version 2 of the License, or
+ ! (at your option) any later version.
+ !
+ ! This program is distributed in the hope that it will be useful,
+ ! but WITHOUT ANY WARRANTY; without even the implied warranty of
+ ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ ! GNU General Public License for more details.
+ !
+ ! You should have received a copy of the GNU General Public License along
+ ! with this program; if not, write to the Free Software Foundation, Inc.,
+ ! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ !
+ !=====================================================================
+ */
+
+/* ---------------------------------------
+
+// asynchronuous file i/o
+
+ --------------------------------------- */
+
+// --------------------------------------------------------------------------------
+// ------------------------- Non-blocking IO routines -----------------------------
+// --------------------------------------------------------------------------------
+
+// In order to overlap disk I/O with computation, we have defined
+// several routines which can do disk I/O in a non-blocking way. The
+// first launches a thread (using pthreads), which writes to
+// disk using fwrite. The second version, attempts to use the linux
+// system call aio_write, which produces a similar result, but is not
+// able to use buffered I/O resulting in slower performance. Both
+// these routines must copy the memory region they are saving, because
+// it will change as the output is uptdated during the next time
+// step. This could be solved by rotating pointers, but this is not
+// very easy in Fortran (especially without changing the fortran code)
+
+
+#include "config.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdbool.h>
+#include <assert.h>
+
+// debug: outputs traces
+#define DEBUG 0
+#if DEBUG == 1
+#define TRACE(x) printf("%s\n",x);
+#else
+#define TRACE(x) // printf("%s\n",x);
+#endif
+
+#define exit_error(msg) { fputs (msg,stderr); abort(); }
+
+// --------------------------------------------------------------------------------
+
+// Pthreads version -- faster than aio version
+
+// --------------------------------------------------------------------------------
+
+// threading version
+#include <pthread.h>
+
+struct threadInfo{
+ int pid;
+ char* buffer;
+ size_t bytes_to_rw;
+ int it_sub;
+ bool finished;
+ bool started;
+};
+
+// --------------------------------------------------------------------------------
+
+// adjoint sources threads
+
+// --------------------------------------------------------------------------------
+
+struct threadInfo ptDataAdj;
+
+pthread_t adj_io_thread;
+
+// fortran function
+extern void FC_FUNC_(read_adjoint_sources_local,READ_ADJOINT_SOURCES_LOCAL)(char*, const int* , const int*);
+
+// Thread that does actual read in adjoint sources. dummy argument is needed to avoid compiler warning.
+void *fread_adj_thread(void* dummy){
+
+ // note: having it_sub_adj as function argument and using int it_sub = (* (int*) it_sub_adj); did not work...
+ TRACE("fread_adj_thread");
+ // debug info
+ //printf("adjoint thread: nadj_rec_local = %i - it_sub_adj = %i \n",ptDataAdj.pid,ptDataAdj.it_sub);
+
+ // calls fortran function
+ // see file: src/specfem3D/read_adjonit_sources.f90
+ FC_FUNC_(read_adjoint_sources_local,READ_ADJOINT_SOURCES_LOCAL)(ptDataAdj.buffer,&ptDataAdj.pid,&ptDataAdj.it_sub);
+
+ // reading done
+ ptDataAdj.finished = true;
+
+ // exit thread
+ pthread_exit(NULL);
+}
+
+// Waits until thread is finished with I/O
+void wait_adj_io_thread() {
+
+ TRACE("wait_adj_io_thread");
+
+ int rc;
+ // checks if thread still runs
+ if( ptDataAdj.started ) {
+ void* status;
+ rc = pthread_join(adj_io_thread, &status);
+ if( rc ) {
+ printf("ERROR; return code from pthread_join() is %d\n", rc);
+ exit_error("error in wait_adj_io_thread: thread_join failed");
+ }
+ // checks finished flag
+ assert(ptDataAdj.finished == true); // Adjoint thread has completed, but somehow it isn't finished?
+
+ // reset
+ ptDataAdj.started = false;
+ }
+}
+
+// initializes adjoint thread
+void
+FC_FUNC_(prepare_adj_io_thread,CREATE_IO_ADJ_THREAD)(char *buffer, long* length,int* nadj_rec_local){
+
+ TRACE("prepare_adj_io_thread");
+
+ int rc;
+ size_t bytes_to_read = *length;
+
+ // checks if buffer valid
+ assert(buffer != NULL); // "Adjoint thread: associated buffer is invalid"
+ if( bytes_to_read <= 0 ) exit_error("Adjoint thread: associated buffer length is invalid");
+
+ // initializes thread info
+ ptDataAdj.started = false;
+ ptDataAdj.finished = false;
+ ptDataAdj.pid = *nadj_rec_local;
+ ptDataAdj.buffer = buffer;
+ ptDataAdj.bytes_to_rw = bytes_to_read;
+ ptDataAdj.it_sub = 0;
+}
+
+// creates thread for reading adjoint sources
+void
+FC_FUNC_(read_adj_io_thread,CREATE_IO_ADJ_THREAD)(int* it_sub_adj){
+
+ TRACE("read_adj_io_thread");
+ // debug
+ //printf("read_adj_io_thread: it_sub_adj = %i \n",*it_sub_adj);
+
+ int rc;
+
+ // checks if buffer valid
+ assert(ptDataAdj.buffer != NULL); // Adjoint thread: associated buffer is invalid
+
+ // waits until previous thread finishes
+ wait_adj_io_thread();
+
+ // prepares the thread
+ pthread_attr_t attr;
+
+ rc = pthread_attr_init(&attr);
+ if( rc != 0 ) exit_error("Adjoint thread: initialization of thread failed");
+
+ rc = pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
+ if( rc != 0 ) exit_error("Adjoint thread: setting thread state failed");
+
+ // sets new thread info
+ ptDataAdj.started = true;
+ ptDataAdj.finished = false;
+ ptDataAdj.it_sub = *it_sub_adj;
+
+ // create and launch the thread.
+ // note: using it_sub_adj as argument (void*) it_sub_adj did not work...
+ rc = pthread_create(&adj_io_thread,&attr,fread_adj_thread,NULL);
+ if( rc != 0 ) exit_error("Adjoint thread: creating thread failed");
+
+}
+
+// synchronizes with read thread to make sure buffer is read in
+void
+FC_FUNC_(sync_adj_io_thread,SYNC_ADJ_IO_THREAD)(char *adj_sourcearrays) {
+
+ TRACE("sync_adj_io_thread");
+
+ // Wait until previous thread finishes reading
+ wait_adj_io_thread();
+
+ // the buffer is updated
+ // copies the thread buffer to the output
+ memcpy(adj_sourcearrays,ptDataAdj.buffer,ptDataAdj.bytes_to_rw);
+}
+
+
+
+// --------------------------------------------------------------------------------
+
+// absorbing boundaries threads
+
+// --------------------------------------------------------------------------------
+// not used yet...
+/*
+struct threadInfo ptData[ABS_FILEID];
+
+pthread_t io_thread[ABS_FILEID];
+
+void *fwrite_thread(void *fileID);
+void *fread_thread(void *fileID);
+void prepare_thread_io(int *fid);
+void write_abs_ptio(int *fid, char *buffer, int *length, int *index);
+void wait_io_thread(int *fid);
+
+// Thread that does actual fwrite.
+void *fwrite_thread(void *fileID)
+{
+ int fid;
+ fid = (int)fileID;
+
+ fwrite(ptData[fid].buffer, 1, ptData[fid].bytes_to_rw,fp_abs[fid]);
+
+ ptData[fid].finished = true;
+ pthread_exit(NULL);
+}
+
+// Thread that does actual fread.
+void *fread_thread(void *fileID)
+{
+ int fid;
+ fid = (int)fileID;
+
+ fread(ptData[fid].buffer, 1, ptData[fid].bytes_to_rw,fp_abs[fid]);
+
+ ptData[fid].finished = true;
+ pthread_exit(NULL);
+}
+
+
+// Setup thread values including the memory region
+void prepare_thread_io(int *fid) {
+ ptData[*fid].started = false;
+ ptData[*fid].finished = false;
+ ptData[*fid].pid = *fid;
+ ptData[*fid].buffer = NULL;
+}
+
+// Starts a thread to write the desired buffer to disk
+void write_abs_ptio(int *fid, char *buffer, int *length, int *index) {
+
+ // allocates buffer if not done so yet
+ size_t bytes_to_write = *length;
+ if(ptData[*fid].buffer == NULL) {
+ ptData[*fid].buffer = (char*)malloc(bytes_to_write);
+ }
+
+ // Wait until previous thread finishes writing
+ wait_io_thread(fid);
+
+ // prepare the thread for writing.
+ pthread_attr_t attr;
+ pthread_attr_init(&attr);
+ pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
+
+ ptData[*fid].started = true;
+ ptData[*fid].finished = false;
+ ptData[*fid].bytes_to_rw = bytes_to_write;
+
+ // the input buffer will be updated while we write the data to disk
+ // -- Thus we have to make a copy so that the data doesn't change
+ // while we are writing it to disk.
+
+ memcpy(ptData[*fid].buffer,buffer,ptData[*fid].bytes_to_rw);
+
+ // create and launch the thread.
+ int rc = pthread_create(&io_thread[*fid],&attr,fwrite_thread,(void*) *fid);
+
+}
+
+// Starts a thread to read the desired buffer from disk
+void read_abs_ptio(int *fid, char *buffer, int *length, int *index) {
+
+ // allocates buffer if not done so yet
+ size_t bytes_to_read = *length;
+ if(ptData[*fid].buffer == NULL) {
+ ptData[*fid].buffer = (char*)malloc(bytes_to_read);
+ }
+
+ // Wait until previous thread finishes reading
+ wait_io_thread(fid);
+
+ // prepare the thread for reading.
+ pthread_attr_t attr;
+ pthread_attr_init(&attr);
+ pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
+
+ ptData[*fid].started = true;
+ ptData[*fid].finished = false;
+ ptData[*fid].bytes_to_rw = bytes_to_read;
+
+ // create and launch the thread.
+ int rc = pthread_create(&io_thread[*fid],&attr,fread_thread,(void*) *fid);
+
+}
+
+
+// Waits until thread is finished with I/O
+void wait_io_thread(int *fid) {
+ int rc;
+ // checks if thread still runs
+ if( ptData[*fid].started ) {
+ void* status;
+ rc = pthread_join(io_thread[*fid], &status);
+ if (rc) {
+ printf("ERROR; return code from pthread_join() is %d\n", rc);
+ exit(-1);
+ }
+ // checks finished flag
+ assert(ptData[*fid].finished == true && "Thread has completed, but somehow it isn't finished?");
+
+ // reset
+ ptData[*fid].started = false;
+ }
+}
+
+// synchronizes with read thread to make sure buffer is read in
+void sync_buffer_io_thread(int *fid, char *buffer, int *length) {
+
+ // checks buffer length
+ size_t bytes_to_read = *length;
+ if(ptData[*fid].buffer == NULL) {
+ assert("Associated thread has no buffer");
+ }
+ assert(ptData[*fid].bytes_to_rw != bytes_to_read && "Associated thread has different buffer length");
+
+ // Wait until previous thread finishes reading
+ wait_io_thread(fid);
+
+ // the input buffer will be updated while we read the data to disk
+ // -- Thus we have to make a copy so that the data doesn't change
+ // while we are reading into it from disk.
+ //
+ // copies the thread buffer to the output
+ memcpy(buffer,ptData[*fid].buffer,ptData[*fid].bytes_to_rw);
+}
+*/
+
+
+
+// --------------------------------------------------------------------------------
+
+// async IO
+
+// --------------------------------------------------------------------------------
+// unused yet...
+/*
+#include <aio.h>
+#include <signal.h>
+#include <errno.h>
+
+// Async_IO
+static int fp_aio[ABS_FILEID];
+struct aiocb fIOw[ABS_FILEID];
+
+// --------------------------------------------------------------------------------
+
+// write routines using AIO
+
+void write_abs_aio(int *fid, char *buffer, int *length, int *index) {
+
+ printf("aio_error return code=%d\n",aio_error(&fIOw[*fid]));
+
+ // wait for previous IO to finish -- aio_error returns 0 when the
+ // job specified by fIOw is finished successfully
+ while(aio_error(&fIOw[*fid])) {}
+
+ // setup async IO for this transfer
+ fIOw[*fid].aio_fildes = fp_aio[*fid];
+ // fIOw[*fid].aio_offset = 0;
+ fIOw[*fid].aio_buf = buffer;
+ size_t write_length = *length;
+ fIOw[*fid].aio_nbytes = write_length;
+ fIOw[*fid].aio_reqprio = 0;
+ fIOw[*fid].aio_sigevent.sigev_notify = SIGEV_NONE;
+ printf("aio_writing %ld bytes\n",write_length);
+ aio_write(&fIOw[*fid]);
+
+}
+
+void close_file_abs_aio(int * fid) {
+
+ // wait for previous IO to finish -- aio_error returns 0 when the
+ // job specified by fIOw is finished successfully
+ while(aio_error(&fIOw[*fid])) {}
+
+ // closes file
+
+ close(fp_aio[*fid]);
+
+}
+*/
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/finalize_simulation.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/finalize_simulation.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -232,6 +232,7 @@
deallocate(iadj_vec)
if(nadj_rec_local > 0) then
deallocate(adj_sourcearrays)
+ if( IO_ASYNC_COPY .and. NSTEP_SUB_ADJ > 1 ) deallocate(buffer_sourcearrays)
deallocate(iadjsrc,iadjsrc_len)
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -346,7 +346,7 @@
z_target(irec) = z_target_rec
! would write out desired target locations of receivers
- !if (myrank == 0) write(IOVTK,*) sngl(x_target(irec)), sngl(y_target(irec)), sngl(z_target(irec))
+ !if (myrank == 0) write(IOUT_VTK,*) sngl(x_target(irec)), sngl(y_target(irec)), sngl(z_target(irec))
! flag to check that we located at least one target element
located_target = .false.
@@ -711,7 +711,7 @@
if(myrank == 0) then
! appends receiver locations to sr.vtk file
- open(IOVTK,file=trim(OUTPUT_FILES)//'/sr_tmp.vtk', &
+ open(IOUT_VTK,file=trim(OUTPUT_FILES)//'/sr_tmp.vtk', &
position='append',status='old',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'Error opening and appending receivers to file sr_tmp.vtk')
@@ -762,13 +762,13 @@
epidist_found(nrec_found) = epidist(irec)
! writes out actual receiver location to vtk file
- write(IOVTK,'(3e18.6)') sngl(x_found(irec)), sngl(y_found(irec)), sngl(z_found(irec))
+ write(IOUT_VTK,'(3e18.6)') sngl(x_found(irec)), sngl(y_found(irec)), sngl(z_found(irec))
endif
enddo
! finishes sr_tmp.vtk file
- write(IOVTK,*)
- close(IOVTK)
+ write(IOUT_VTK,*)
+ close(IOUT_VTK)
! compute maximal distance for all the receivers
final_distance_max = maxval(final_distance(:))
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -179,7 +179,7 @@
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
- open(IOVTK,file=trim(OUTPUT_FILES)//'/sr_tmp.vtk', &
+ open(IOUT_VTK,file=trim(OUTPUT_FILES)//'/sr_tmp.vtk', &
position='append',status='old',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'Error opening and appending sources to file sr_tmp.vtk')
endif
@@ -669,7 +669,7 @@
write(IMAIN,*) ' time shift: ',tshift_cmt(isource),' seconds'
! writes out actual source position to vtk file
- write(IOVTK,'(3e18.6)') sngl(x_found_source(isource_in_this_subset)), &
+ write(IOUT_VTK,'(3e18.6)') sngl(x_found_source(isource_in_this_subset)), &
sngl(y_found_source(isource_in_this_subset)), &
sngl(z_found_source(isource_in_this_subset))
@@ -747,7 +747,7 @@
call flush_IMAIN()
! closing sr_tmp.vtk
- close(IOVTK)
+ close(IOUT_VTK)
endif
call synchronize_all()
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_adjoint_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_adjoint_sources.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_adjoint_sources.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -0,0 +1,265 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and CNRS / INRIA / University of Pau, France
+! (c) Princeton University and CNRS / INRIA / University of Pau
+! August 2013
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine read_adjoint_sources()
+
+! reads in adjoint source files
+!
+! note: depending on the number of adjoint source files and cluster filesystem,
+! this can take a very long time and slow down adjoint simulations considerably.
+!
+! when flag IO_ASYNC_COPY is set, a new parallel thread is created to read in adjoint source files,
+! while the main processes can continue computing. we thus overlap file i/o with computations.
+! in case the heavy computations take long enough, this overlap should perfectly hide the i/o latency.
+!
+! the costs are then for the additional memory, i.e. buffer array (buffer_sourcearrays), and the
+! copying operation of data from the buffer array to the actual adj_sourcearrays array.
+
+ use specfem_par
+
+ implicit none
+
+ ! local parameters
+ integer :: it_sub_adj
+ ! debug timing
+ !double precision, external :: wtime
+ !double precision :: tstart
+
+
+ ! we already checked that nadj_rec_local > 0 before calling this routine
+
+ ! debug timing
+ !tstart = wtime()
+
+ ! determines chunk_number
+ it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) )
+
+ ! debug
+ !print*,'read adjoint sources: it_sub_adj = ',it_sub_adj
+
+ ! asynchronuously reads in adjoint source files
+ if( IO_ASYNC_COPY .and. NSTEP_SUB_ADJ > 1 ) then
+ ! handles file input/output thread
+ if( it == it_begin ) then
+ ! creates new io thread for reading in sources
+ call read_adj_io_thread(it_sub_adj)
+
+ ! first chunk of adjoint sources must ready at begining, so we wait.
+ ! waits for previous read to finish and
+ ! copy over buffered data into tmp_sourcearray
+ call sync_adj_io_thread(adj_sourcearrays)
+
+ else
+ ! waits for previous read to finish and
+ ! copy over buffered data into tmp_sourcearray
+ call sync_adj_io_thread(adj_sourcearrays)
+ endif
+
+ ! checks if next chunk necessary
+ if( it_sub_adj < NSTEP_SUB_ADJ ) then
+ ! starts thread to read in next junk
+ it_sub_adj = it_sub_adj + 1
+
+ ! creates new io thread for reading in sources
+ call read_adj_io_thread(it_sub_adj)
+ endif
+
+ else
+ ! synchronuous read routine
+
+ ! reads in local adjoint sources
+ call read_adjoint_sources_local(adj_sourcearrays,nadj_rec_local,it_sub_adj)
+
+ endif
+
+ ! debug timing
+ !print*,'read adjoint sources: elapsed time = ',wtime() - tstart
+ !print*
+
+ end subroutine read_adjoint_sources
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine read_adjoint_sources_local(sourcearrays,nadj_rec_local,it_sub_adj)
+
+! reads in local adjoint source files
+
+ use specfem_par,only: myrank,NPROCTOT_VAL, &
+ nrec,islice_selected_rec,station_name,network_name, &
+ xi_receiver,eta_receiver,gamma_receiver,nu,xigll,yigll,zigll, &
+ iadjsrc_len,iadjsrc,NSTEP_SUB_ADJ, &
+ DT,CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ,NTSTEP_BETWEEN_READ_ADJSRC
+
+ implicit none
+
+ integer,intent(in) :: nadj_rec_local
+ integer,intent(in) :: it_sub_adj
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ,nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC) :: sourcearrays
+
+ ! local parameters
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: tmp_sourcearray
+ integer :: irec,irec_local,itime,ier
+ character(len=150) :: adj_source_file
+
+ ! debug
+ !print*,'reading adjoint sources local:',myrank,' - chunk ',it_sub_adj,'out of ',NSTEP_SUB_ADJ, &
+ ! ' for local adjoint sources = ',nadj_rec_local
+
+ ! checks chunk number
+ if( it_sub_adj < 1 .or. it_sub_adj > NSTEP_SUB_ADJ ) then
+ print*,'error reading adjoint sources: chunk number ',it_sub_adj,'is invalid'
+ call exit_MPI(myrank,'error reading adjoint sources with invalid chunk number')
+ endif
+
+ ! allocates temporary source array
+ allocate(tmp_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NTSTEP_BETWEEN_READ_ADJSRC),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating array tmp_sourcearray')
+
+ ! initializes
+ tmp_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
+
+ ! counter
+ irec_local = 0
+
+ ! loops over all adjoint sources
+ do irec = 1, nrec
+ ! checks that the source slice number is okay
+ if(islice_selected_rec(irec) < 0 .or. islice_selected_rec(irec) > NPROCTOT_VAL-1) then
+ print*,'error rank ',myrank,': adjoint source slice index ',islice_selected_rec(irec),&
+ ' is out of bounds ',NPROCTOT_VAL-1
+ call exit_MPI(myrank,'error adjoint source has wrong source slice number in adjoint simulation')
+ endif
+
+ ! compute source arrays for adjoint sources within this rank's slice
+ if( myrank == islice_selected_rec(irec) ) then
+ ! increases counter
+ irec_local = irec_local + 1
+
+ ! adjoint source file name **sta**.**net**
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+
+ ! reads in **sta**.**net**.**.adj files
+ call compute_arrays_source_adjoint(myrank,adj_source_file, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ nu(:,:,irec),tmp_sourcearray, &
+ xigll,yigll,zigll, &
+ iadjsrc_len(it_sub_adj),iadjsrc,it_sub_adj, &
+ NSTEP_SUB_ADJ,NTSTEP_BETWEEN_READ_ADJSRC,DT)
+
+ ! stores source array
+ ! note: the adj_sourcearrays has a time stepping from 1 to NTSTEP_BETWEEN_READ_ADJSRC
+ ! this gets overwritten every time a new block/chunk is read in
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ sourcearrays(:,:,:,:,irec_local,itime) = tmp_sourcearray(:,:,:,:,itime)
+ enddo
+
+ endif
+ enddo
+
+ ! checks that number of read sources is valid
+ if(irec_local /= nadj_rec_local) then
+ call exit_MPI(myrank,'irec_local /= nadj_rec_local in adjoint simulation')
+ endif
+
+ ! frees temporary array
+ deallocate(tmp_sourcearray)
+
+ end subroutine read_adjoint_sources_local
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine check_adjoint_sources(irec,nadj_files_found)
+
+ use specfem_par
+
+ implicit none
+
+ integer,intent(in) :: irec
+ ! counter
+ integer,intent(inout) :: nadj_files_found
+
+ ! local parameters
+ double precision :: junk
+ integer :: icomp,itime
+ integer :: ier
+ character(len=256) :: filename,adj_source_file
+ character(len=3),dimension(NDIM) :: comp
+ character(len=2) :: bic
+
+ ! root file name
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+
+ ! bandwidth code
+ ! by Ebru
+ call band_instrument_code(DT,bic)
+ comp(1) = bic(1:2)//'N'
+ comp(2) = bic(1:2)//'E'
+ comp(3) = bic(1:2)//'Z'
+
+ ! loops over file components E/N/Z
+ do icomp = 1,NDIM
+
+ ! opens adjoint source file for this component
+ filename = 'SEM/'//trim(adj_source_file) // '.'// comp(icomp) // '.adj'
+ open(unit=IIN,file=trim(filename),status='old',action='read',iostat=ier)
+
+ ! checks if file opens/exists
+ if( ier /= 0 ) then
+ ! adjoint source file not found
+ ! stops simulation
+ call exit_MPI(myrank,&
+ 'file '//trim(filename)//' not found, please check with your STATIONS_ADJOINT file')
+ endif
+
+ ! checks length of file
+ itime = 0
+ do while(ier == 0)
+ read(IIN,*,iostat=ier) junk,junk
+ if( ier == 0 ) itime = itime + 1
+ enddo
+
+ ! checks length
+ if( itime /= NSTEP) then
+ print*,'adjoint source error: ',trim(filename),' has length',itime,' but should be',NSTEP
+ call exit_MPI(myrank,&
+ 'file '//trim(filename)//' length is wrong, please check your adjoint sources and your simulation duration')
+ endif
+
+ ! updates counter for found files
+ nadj_files_found = nadj_files_found + 1
+
+ ! closes file
+ close(IIN)
+
+ enddo
+
+ end subroutine check_adjoint_sources
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk 2014-01-24 16:03:04 UTC (rev 22994)
@@ -40,6 +40,7 @@
$O/compute_arrays_source.solver.o \
$O/convert_time.solver.o \
$O/define_derivation_matrices.solver.o \
+ $O/file_io_threads.cc.o \
$O/get_backazimuth.solver.o \
$O/get_cmt.solver.o \
$O/get_event_info.solver.o \
@@ -81,6 +82,7 @@
$O/multiply_arrays_source.solverstatic.o \
$O/noise_tomography.solverstatic.o \
$O/prepare_timerun.solverstatic.o \
+ $O/read_adjoint_sources.solverstatic.o \
$O/read_arrays_solver.solverstatic.o \
$O/read_forward_arrays.solverstatic.o \
$O/read_mesh_databases.solverstatic.o \
@@ -335,6 +337,9 @@
$O/%.solver.o: $S/%.F90 $O/shared_par.shared_module.o
${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $<
+$O/%.cc.o: $S/%.c ${SETUP}/config.h
+ ${CC} -c $(CPPFLAGS) $(CFLAGS) -o $@ $<
+
###
### CUDA 5 only
###
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -126,16 +126,16 @@
if(myrank == 0) then
! write source and receiver VTK files for Paraview
filename = trim(OUTPUT_FILES)//'/sr_tmp.vtk'
- open(IOVTK,file=trim(filename),status='unknown',iostat=ier)
+ open(IOUT_VTK,file=trim(filename),status='unknown',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error opening temporary file sr_temp.vtk')
- write(IOVTK,'(a)') '# vtk DataFile Version 2.0'
- write(IOVTK,'(a)') 'Source and Receiver VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 2.0'
+ write(IOUT_VTK,'(a)') 'Source and Receiver VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
! LQY -- won't be able to know NSOURCES+nrec at this point...
- write(IOVTK, '(a,i6,a)') 'POINTS ', NSOURCES, ' float'
+ write(IOUT_VTK, '(a,i6,a)') 'POINTS ', NSOURCES, ' float'
! closing file, rest of informations will be appended later on
- close(IOVTK)
+ close(IOUT_VTK)
endif
! locate sources in the mesh
@@ -331,12 +331,8 @@
implicit none
! local parameters
- double precision :: junk
integer :: irec,isource,nrec_tot_found
- integer :: icomp,itime,nadj_files_found,nadj_files_found_tot
- character(len=3),dimension(NDIM) :: comp
- character(len=256) :: filename,adj_source_file
- character(len=2) :: bic
+ integer :: nadj_files_found,nadj_files_found_tot
integer :: ier
integer,dimension(:),allocatable :: tmp_rec_local_all
integer :: maxrec,maxproc(1)
@@ -407,12 +403,6 @@
! counts receivers for adjoint simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
- ! by Ebru
- call band_instrument_code(DT,bic)
- comp(1) = bic(1:2)//'N'
- comp(2) = bic(1:2)//'E'
- comp(3) = bic(1:2)//'Z'
-
! temporary counter to check if any files are found at all
nadj_files_found = 0
do irec = 1,nrec
@@ -425,37 +415,7 @@
nadj_rec_local = nadj_rec_local + 1
! checks **sta**.**net**.**MX**.adj files for correct number of time steps
- adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
- do icomp = 1,NDIM
-
- ! opens adjoint source file for this component
- filename = 'SEM/'//trim(adj_source_file) // '.'// comp(icomp) // '.adj'
- open(unit=IIN,file=trim(filename),status='old',action='read',iostat=ier)
-
- if( ier == 0 ) then
- ! checks length of file
- itime = 0
- do while(ier == 0)
- read(IIN,*,iostat=ier) junk,junk
- if( ier == 0 ) itime = itime + 1
- enddo
- ! checks length
- if( itime /= NSTEP) then
- print*,'adjoint source error: ',trim(filename),' has length',itime,' but should be',NSTEP
- call exit_MPI(myrank,&
- 'file '//trim(filename)//' has wrong length, please check with your simulation duration')
- endif
-
- ! updates counter for found files
- nadj_files_found = nadj_files_found + 1
- else
- ! adjoint source file not found
- ! stops simulation
- call exit_MPI(myrank,&
- 'file '//trim(filename)//' not found, please check with your STATIONS_ADJOINT file')
- endif
- close(IIN)
- enddo
+ call check_adjoint_sources(irec,nadj_files_found)
endif
enddo
@@ -481,6 +441,7 @@
else
write(IMAIN,*) 'this total is okay'
endif
+ write(IMAIN,*)
call flush_IMAIN()
endif
@@ -502,6 +463,24 @@
if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary array tmp_rec_local_all')
endif
+ ! user output infos
+ ! sources
+ if( SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3 ) then
+ ! user output
+ if( myrank == 0 ) then
+ ! note: all process allocate the full sourcearrays array
+ ! sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,NSOURCES)
+ size = dble(NSOURCES) * dble(NDIM * NGLLX * NGLLY * NGLLZ * CUSTOM_REAL / 1024. / 1024. )
+ ! outputs info
+ write(IMAIN,*) 'source arrays:'
+ write(IMAIN,*) ' number of sources is ',NSOURCES
+ write(IMAIN,*) ' size of source array = ', sngl(size),'MB'
+ write(IMAIN,*) ' = ', sngl(size/1024.d0),'GB'
+ write(IMAIN,*)
+ call flush_IMAIN()
+ endif
+ endif
+
! seismograms
! gather from slaves on master
tmp_rec_local_all(:) = 0
@@ -523,12 +502,11 @@
size = dble(maxrec) * dble(NDIM * NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * CUSTOM_REAL / 1024. / 1024. )
endif
! outputs info
- write(IMAIN,*)
write(IMAIN,*) 'seismograms:'
write(IMAIN,*) ' writing out seismograms at every NTSTEP_BETWEEN_OUTPUT_SEISMOS = ',NTSTEP_BETWEEN_OUTPUT_SEISMOS
write(IMAIN,*) ' maximum number of local receivers is ',maxrec,' in slice ',maxproc(1)
- write(IMAIN,*) ' size of maximum seismogram array = ', sngl(size),'MB'
- write(IMAIN,*) ' = ', sngl(size/1024.d0),'GB'
+ write(IMAIN,*) ' size of maximum seismogram array = ', sngl(size),'MB'
+ write(IMAIN,*) ' = ', sngl(size/1024.d0),'GB'
write(IMAIN,*)
call flush_IMAIN()
endif
@@ -555,10 +533,15 @@
! adj_sourcearrays size in MB
! adj_sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC)
size = dble(maxrec) * dble(NDIM * NGLLX * NGLLY * NGLLZ * NTSTEP_BETWEEN_READ_ADJSRC * CUSTOM_REAL / 1024. / 1024. )
-
+ ! note: in case IO_ASYNC_COPY is set, and depending of NSTEP_SUB_ADJ,
+ ! this memory requirement might double.
+ ! at this point, NSTEP_SUB_ADJ is not set yet...
! outputs info
write(IMAIN,*) 'adjoint source arrays:'
write(IMAIN,*) ' reading adjoint sources at every NTSTEP_BETWEEN_READ_ADJSRC = ',NTSTEP_BETWEEN_READ_ADJSRC
+ if( IO_ASYNC_COPY ) then
+ write(IMAIN,*) ' using asynchronuous buffer for file i/o of adjoint sources'
+ endif
write(IMAIN,*) ' maximum number of local adjoint sources is ',maxrec,' in slice ',maxproc(1)
write(IMAIN,*) ' size of maximum adjoint source array = ', sngl(size),'MB'
write(IMAIN,*) ' = ', sngl(size/1024.d0),'GB'
@@ -569,7 +552,6 @@
deallocate(tmp_rec_local_all)
-
end subroutine setup_receivers
!
@@ -627,6 +609,7 @@
! local parameters
integer :: ier
+ integer(kind=8) :: arraysize
! allocates source arrays
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
@@ -650,7 +633,7 @@
iadj_vec(it) = NSTEP-it+1 ! default is for reversing entire record, e.g. 3000,2999,..,1
enddo
- ! number of adjoint source blocks to read in
+ ! total number of adjoint source blocks to read in
NSTEP_SUB_ADJ = ceiling( dble(NSTEP)/dble(NTSTEP_BETWEEN_READ_ADJSRC) )
if(nadj_rec_local > 0) then
@@ -660,6 +643,24 @@
if( ier /= 0 ) call exit_MPI(myrank,'error allocating adjoint sourcearrays')
adj_sourcearrays(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ ! additional buffer for asynchronuous file i/o
+ if( IO_ASYNC_COPY .and. NSTEP_SUB_ADJ > 1 ) then
+ ! allocates read buffer
+ allocate(buffer_sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC), &
+ stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating array buffer_sourcearrays')
+
+ ! array size in bytes (note: the multiplication is split into two line to avoid integer-overflow)
+ arraysize = NDIM * NGLLX * NGLLY * NGLLZ * CUSTOM_REAL
+ arraysize = arraysize * nadj_rec_local * NTSTEP_BETWEEN_READ_ADJSRC
+
+ ! debug
+ !print*,'buffer_sourcearrays: size = ',arraysize,' Bytes = ',arraysize/1024./1024.,'MB'
+
+ ! initializes io thread
+ call prepare_adj_io_thread(buffer_sourcearrays,arraysize,nadj_rec_local)
+ endif
+
! allocate indexing arrays
allocate(iadjsrc(NSTEP_SUB_ADJ,2), &
iadjsrc_len(NSTEP_SUB_ADJ),stat=ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -166,12 +166,17 @@
! Lagrange interpolators at receivers
double precision, dimension(:,:), allocatable :: hxir_store,hetar_store,hgammar_store
- !ADJOINT
+ ! ADJOINT sources
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:,:), allocatable :: adj_sourcearrays
+ ! asynchronuous read buffer when IO_ASYNC_COPY is set
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:,:), allocatable :: buffer_sourcearrays
+
integer :: nrec_simulation, nadj_rec_local
integer :: NSTEP_SUB_ADJ ! to read input in chunks
+
integer, dimension(:,:), allocatable :: iadjsrc ! to read input in chunks
integer, dimension(:), allocatable :: iadjsrc_len,iadj_vec
+
! source frechet derivatives
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: moment_der
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: sloc_der
@@ -179,6 +184,7 @@
double precision, dimension(:,:), allocatable :: hpxir_store,hpetar_store,hpgammar_store
integer :: nadj_hprec_local
+
!-----------------------------------------------------------------
! seismograms
!-----------------------------------------------------------------
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -32,6 +32,8 @@
implicit none
+ ! note: this routine gets called if( nrec_local > 0 .or. ( WRITE_SEISMOGRAMS_BY_MASTER .and. myrank == 0 ) )
+
! update position in seismograms
seismo_current = seismo_current + 1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Visualization/VTK_ParaView/create_slice_VTK.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Visualization/VTK_ParaView/create_slice_VTK.f90 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Visualization/VTK_ParaView/create_slice_VTK.f90 2014-01-24 16:03:04 UTC (rev 22994)
@@ -268,86 +268,86 @@
write(IMAIN,*) ' vtk file: '
write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
- open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
- write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
- write(IOVTK,'(a)') 'material model VTK file'
- write(IOVTK,'(a)') 'ASCII'
- write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ open(IOUT_VTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOUT_VTK,'(a)') 'material model VTK file'
+ write(IOUT_VTK,'(a)') 'ASCII'
+ write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
! writes out all points for each element, not just global ones
- write(IOVTK, '(a,i12,a)') 'POINTS ', nspec*8, ' float'
+ write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nspec*8, ' float'
do ispec=1,nspec
i = ibool(1,1,1,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(NGLLX,1,1,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(NGLLX,NGLLY,1,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(1,NGLLY,1,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(1,1,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(NGLLX,1,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(NGLLX,NGLLY,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
i = ibool(1,NGLLY,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! note: indices for vtk start at 0
- write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
do ispec=1,nspec
- write(IOVTK,'(9i12)') 8,(ispec-1)*8,(ispec-1)*8+1,(ispec-1)*8+2,(ispec-1)*8+3,&
+ write(IOUT_VTK,'(9i12)') 8,(ispec-1)*8,(ispec-1)*8+1,(ispec-1)*8+2,(ispec-1)*8+3,&
(ispec-1)*8+4,(ispec-1)*8+5,(ispec-1)*8+6,(ispec-1)*8+7
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
! type: hexahedrons
- write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
- write(IOVTK,'(6i12)') (12,ispec=1,nspec)
- write(IOVTK,*) ""
+ write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOUT_VTK,'(6i12)') (12,ispec=1,nspec)
+ write(IOUT_VTK,*) ""
! writes out gll-data (velocity) for each element point
- write(IOVTK,'(a,i12)') "POINT_DATA ",nspec*8
- write(IOVTK,'(a)') "SCALARS gll_data float"
- write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ write(IOUT_VTK,'(a,i12)') "POINT_DATA ",nspec*8
+ write(IOUT_VTK,'(a)') "SCALARS gll_data float"
+ write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
do ispec = 1,nspec
i = ibool(1,1,1,ispec)
- write(IOVTK,'(3e18.6)') gll_data(1,1,1,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(1,1,1,ispec)
i = ibool(NGLLX,1,1,ispec)
- write(IOVTK,'(3e18.6)') gll_data(NGLLX,1,1,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(NGLLX,1,1,ispec)
i = ibool(NGLLX,NGLLY,1,ispec)
- write(IOVTK,'(3e18.6)') gll_data(NGLLX,NGLLY,1,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(NGLLX,NGLLY,1,ispec)
i = ibool(1,NGLLY,1,ispec)
- write(IOVTK,'(3e18.6)') gll_data(1,NGLLY,1,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(1,NGLLY,1,ispec)
i = ibool(1,1,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') gll_data(1,1,NGLLZ,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(1,1,NGLLZ,ispec)
i = ibool(NGLLX,1,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') gll_data(NGLLX,1,NGLLZ,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(NGLLX,1,NGLLZ,ispec)
i = ibool(NGLLX,NGLLY,NGLLZ,ispec)
- write(IOVTK,'(3e18.6)') gll_data(NGLLX,NGLLY,NGLLZ,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(NGLLX,NGLLY,NGLLZ,ispec)
i = ibool(1,NGLLY,NGLLZ,ispec)-1
- write(IOVTK,'(3e18.6)') gll_data(1,NGLLY,NGLLZ,ispec)
+ write(IOUT_VTK,'(3e18.6)') gll_data(1,NGLLY,NGLLZ,ispec)
enddo
- write(IOVTK,*) ""
+ write(IOUT_VTK,*) ""
- close(IOVTK)
+ close(IOUT_VTK)
end subroutine write_VTK_data_gll_cr
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v36.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v36.h 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v36.h 2014-01-24 16:03:04 UTC (rev 22994)
@@ -95,7 +95,7 @@
! uncomment this to write messages to the screen (slows down the code)
! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
! I/O unit for source and receiver vtk file
- integer, parameter :: IOVTK = 98
+ integer, parameter :: IOUT_VTK = 98
! R_EARTH is the radius of the bottom of the oceans (radius of Earth in m)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v40.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v40.h 2014-01-21 15:03:43 UTC (rev 22993)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/doubling_brick/constants_modified_v40.h 2014-01-24 16:03:04 UTC (rev 22994)
@@ -95,7 +95,7 @@
! uncomment this to write messages to the screen (slows down the code)
! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
! I/O unit for source and receiver vtk file
- integer, parameter :: IOVTK = 98
+ integer, parameter :: IOUT_VTK = 98
! R_EARTH is the radius of the bottom of the oceans (radius of Earth in m)
More information about the CIG-COMMITS
mailing list