[cig-commits] r20319 - in seismo/2D/SPECFEM2D/trunk/src: meshfem2D specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Jun 5 09:17:12 PDT 2012
Author: dkomati1
Date: 2012-06-05 09:17:11 -0700 (Tue, 05 Jun 2012)
New Revision: 20319
Modified:
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90
Log:
subsamp_seismos now implemented
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90 2012-06-05 14:01:31 UTC (rev 20318)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_parameter_file.F90 2012-06-05 16:17:11 UTC (rev 20319)
@@ -280,6 +280,7 @@
call read_value_integer_p(subsamp_seismos, 'solver.subsamp_seismos')
if(err_occurred() /= 0) stop 'error reading parameter 33e in Par_file'
+ if(subsamp_seismos < 1) stop 'error: subsamp_seismos must be >= 1'
call read_value_logical_p(generate_STATIONS, 'solver.generate_STATIONS')
if(err_occurred() /= 0) stop 'error reading parameter 23 in Par_file'
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2012-06-05 14:01:31 UTC (rev 20318)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2012-06-05 16:17:11 UTC (rev 20319)
@@ -295,13 +295,27 @@
read(IIN,*) NSTEP,deltat
if (myrank == 0 .and. ipass == 1) write(IOUT,703) NSTEP,deltat,NSTEP*deltat
- if( SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. &
- (ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) ) then
- print*, '*************** WARNING ***************'
- print*, 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
- stop
+ if(SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. &
+ (ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART)) then
+ print *, '*************** error ***************'
+ stop 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
endif
+! make sure NSTEP_BETWEEN_OUTPUT_SEISMOS is a multiple of subsamp_seismos
+ if(mod(NSTEP_BETWEEN_OUTPUT_SEISMOS,subsamp_seismos) /= 0) &
+ stop 'error: NSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of subsamp_seismos'
+
+! make sure NSTEP is a multiple of subsamp_seismos
+! if not, increase it a little bit, to the next multiple
+ if(mod(NSTEP,subsamp_seismos) /= 0) then
+ if (myrank == 0) then
+ print *,'NSTEP is not a multiple of subsamp_seismos'
+ NSTEP = (NSTEP/subsamp_seismos + 1)*subsamp_seismos
+ print *,'thus increasing it automatically to the next multiple, which is ',NSTEP
+ print *
+ endif
+ endif
+
! output seismograms at least once at the end of the simulation
NSTEP_BETWEEN_OUTPUT_SEISMOS = min(NSTEP,NSTEP_BETWEEN_OUTPUT_SEISMOS)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-06-05 14:01:31 UTC (rev 20318)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-06-05 16:17:11 UTC (rev 20319)
@@ -2355,9 +2355,9 @@
! allocate seismogram arrays
if(ipass == 1) then
- allocate(sisux(NSTEP_BETWEEN_OUTPUT_SEISMOS,nrecloc))
- allocate(sisuz(NSTEP_BETWEEN_OUTPUT_SEISMOS,nrecloc))
- allocate(siscurl(NSTEP_BETWEEN_OUTPUT_SEISMOS,nrecloc))
+ allocate(sisux(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc))
+ allocate(sisuz(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc))
+ allocate(siscurl(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc))
endif
! check if acoustic receiver is exactly on the free surface because pressure is zero there
@@ -3971,7 +3971,7 @@
endif
#ifdef USE_MPI
- if(output_energy) stop 'energy calculation currently serial only, should add an MPI_REDUCE in parallel'
+ if(output_energy) stop 'error: energy calculation currently serial only, should add an MPI_REDUCE in parallel'
#endif
!---- create a Gnuplot script to display the energy curve in log scale
@@ -4340,9 +4340,6 @@
do it = 1,NSTEP
-! update position in seismograms
- seismo_current = seismo_current + 1
-
! compute current time
time = (it-1)*deltat
@@ -7067,7 +7064,6 @@
!---- display time step and max of norm of displacement
if(mod(it,NSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
-
call check_stability(myrank,time,it,NSTEP,NOISE_TOMOGRAPHY, &
nglob_acoustic,nglob_elastic,nglob_poroelastic, &
any_elastic_glob,any_elastic,displ_elastic, &
@@ -7080,10 +7076,14 @@
! add a barrier if we generate traces of the run for analysis with "ParaVer"
if(GENERATE_PARAVER_TRACES) call MPI_BARRIER(MPI_COMM_WORLD,ier)
#endif
-
endif
!---- loop on all the receivers to compute and store the seismograms
+ if(mod(it-1,subsamp_seismos) == 0) then
+
+! update position in seismograms
+ seismo_current = seismo_current + 1
+
do irecloc = 1,nrecloc
irec = recloc(irecloc)
@@ -7223,6 +7223,10 @@
enddo
enddo
+ ! check for edge effects
+ if(seismo_current < 1 .or. seismo_current > NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos) &
+ stop 'error: seismo_current out of bounds in recording of seismograms'
+
! rotate seismogram components if needed, except if recording pressure, which is a scalar
if(seismotype /= 4 .and. seismotype /= 6) then
if(p_sv) then
@@ -7240,6 +7244,7 @@
enddo
+ endif
!----- writing the kernels
! kernels output
@@ -8021,7 +8026,7 @@
nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current,p_sv, &
st_zval,x_source(1),z_source(1),SU_FORMAT,save_ASCII_seismograms, &
- save_binary_seismograms_single,save_binary_seismograms_double)
+ save_binary_seismograms_single,save_binary_seismograms_double,subsamp_seismos)
seismo_offset = seismo_offset + seismo_current
seismo_current = 0
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90 2012-06-05 14:01:31 UTC (rev 20318)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90 2012-06-05 16:17:11 UTC (rev 20319)
@@ -48,7 +48,7 @@
NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current,p_sv, &
st_zval,x_source,z_source,SU_FORMAT,save_ASCII_seismograms, &
- save_binary_seismograms_single,save_binary_seismograms_double)
+ save_binary_seismograms_single,save_binary_seismograms_double,subsamp_seismos)
implicit none
@@ -57,7 +57,7 @@
include "mpif.h"
#endif
- integer :: nrec,NSTEP,seismotype
+ integer :: nrec,NSTEP,seismotype,subsamp_seismos
integer :: NSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current
double precision :: t0,deltat
@@ -69,7 +69,7 @@
integer, intent(in) :: nrecloc,myrank
integer, dimension(nrec),intent(in) :: which_proc_receiver
- double precision, dimension(NSTEP_BETWEEN_OUTPUT_SEISMOS,nrecloc), intent(in) :: sisux,sisuz,siscurl
+ double precision, dimension(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc), intent(in) :: sisux,sisuz,siscurl
double precision :: st_xval(nrec)
@@ -136,7 +136,7 @@
number_of_components = NDIM
endif
- allocate(buffer_binary(NSTEP_BETWEEN_OUTPUT_SEISMOS,number_of_components))
+ allocate(buffer_binary(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,number_of_components))
if (save_binary_seismograms .and. myrank == 0 .and. seismo_offset == 0) then
@@ -226,16 +226,16 @@
#ifdef USE_MPI
else
- call MPI_RECV(buffer_binary(1,1),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,&
+ call MPI_RECV(buffer_binary(1,1),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
if ( number_of_components == 2 ) then
- call MPI_RECV(buffer_binary(1,2),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,&
+ call MPI_RECV(buffer_binary(1,2),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
endif
if ( number_of_components == 3 ) then
- call MPI_RECV(buffer_binary(1,2),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,&
+ call MPI_RECV(buffer_binary(1,2),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
- call MPI_RECV(buffer_binary(1,3),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,&
+ call MPI_RECV(buffer_binary(1,3),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
endif
#endif
@@ -378,14 +378,17 @@
#ifdef USE_MPI
else
- if ( which_proc_receiver(irec) == myrank ) then
+ if (which_proc_receiver(irec) == myrank) then
irecloc = irecloc + 1
- call MPI_SEND(sisux(1,irecloc),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
- if ( number_of_components >= 2 ) then
- call MPI_SEND(sisuz(1,irecloc),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ call MPI_SEND(sisux(1,irecloc),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,0,irec, &
+ MPI_COMM_WORLD,ierror)
+ if (number_of_components >= 2) then
+ call MPI_SEND(sisuz(1,irecloc),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,0,irec, &
+ MPI_COMM_WORLD,ierror)
endif
- if ( number_of_components == 3 ) then
- call MPI_SEND(siscurl(1,irecloc),NSTEP_BETWEEN_OUTPUT_SEISMOS,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ if (number_of_components == 3) then
+ call MPI_SEND(siscurl(1,irecloc),NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,MPI_DOUBLE_PRECISION,0,irec, &
+ MPI_COMM_WORLD,ierror)
endif
endif
#endif
More information about the CIG-COMMITS
mailing list