[cig-commits] [commit] devel: updates get_cmt() and get_force(); only master process reads in source files (49ec9de)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Sep 25 03:49:43 PDT 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/69d1a6fe5d756c2d8d122dd4521d96ddd99b5671...63b33a6cfd79d58ef4ebcf2ca10eb8317b416a7d
>---------------------------------------------------------------
commit 49ec9def6a70f69c9eaa77845d5a0c388e49682c
Author: daniel peter <peterda at ethz.ch>
Date: Mon Sep 22 11:38:17 2014 +0200
updates get_cmt() and get_force(); only master process reads in source files
>---------------------------------------------------------------
49ec9def6a70f69c9eaa77845d5a0c388e49682c
src/generate_databases/rules.mk | 2 -
src/shared/get_cmt.f90 | 362 +++++++++++++++++++++++++++---
src/shared/get_force.f90 | 47 ++--
src/shared/parallel.f90 | 19 ++
src/shared/serial.f90 | 16 ++
src/specfem3D/locate_source.f90 | 53 ++++-
src/specfem3D/prepare_timerun.F90 | 157 +++++++++++++
src/specfem3D/setup_sources_receivers.f90 | 1 -
8 files changed, 589 insertions(+), 68 deletions(-)
diff --git a/src/generate_databases/rules.mk b/src/generate_databases/rules.mk
index 818a844..9b7918e 100644
--- a/src/generate_databases/rules.mk
+++ b/src/generate_databases/rules.mk
@@ -95,9 +95,7 @@ generate_databases_SHARED_OBJECTS = \
$O/detect_surface.shared.o \
$O/exit_mpi.shared.o \
$O/get_attenuation_model.shared.o \
- $O/get_cmt.shared.o \
$O/get_element_face.shared.o \
- $O/get_force.shared.o \
$O/get_global.shared.o \
$O/get_jacobian_boundaries.shared.o \
$O/get_shape2D.shared.o \
diff --git a/src/shared/get_cmt.f90 b/src/shared/get_cmt.f90
index c2f0261..5b6dfe0 100644
--- a/src/shared/get_cmt.f90
+++ b/src/shared/get_cmt.f90
@@ -28,7 +28,7 @@
subroutine get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor,&
DT,NSOURCES,min_tshift_cmt_original)
- use constants
+ use constants,only: IIN,IN_DATA_FILES_PATH,MAX_STRING_LEN,NUMBER_OF_SIMULTANEOUS_RUNS,mygroup
implicit none
@@ -42,12 +42,13 @@
double precision, dimension(NSOURCES), intent(out) :: tshift_cmt,hdur,lat,long,depth
double precision, dimension(6,NSOURCES), intent(out) :: moment_tensor
-!--- local variables below
-
- integer mo,da,julian_day,isource
- double precision t_shift(NSOURCES)
- character(len=5) datasource
- character(len=MAX_STRING_LEN) :: string, CMTSOLUTION,path_to_add
+ ! local variables below
+ integer :: mo,da,julian_day,isource
+ integer :: i,itype,istart,iend,ier
+ double precision :: t_shift(NSOURCES)
+ !character(len=5) :: datasource
+ character(len=256) :: string
+ character(len=MAX_STRING_LEN) :: CMTSOLUTION,path_to_add
! initializes
lat(:) = 0.d0
@@ -57,11 +58,6 @@
tshift_cmt(:) = 0.d0
hdur(:) = 0.d0
moment_tensor(:,:) = 0.d0
- yr = 0
- jda = 0
- ho = 0
- mi = 0
- sec = 0.d0
!
!---- read hypocenter info
@@ -75,76 +71,251 @@
CMTSOLUTION = path_to_add(1:len_trim(path_to_add))//CMTSOLUTION(1:len_trim(CMTSOLUTION))
endif
- open(unit=1,file=CMTSOLUTION,status='old',action='read')
+ open(unit = IIN,file=trim(CMTSOLUTION),status='old',action='read',iostat=ier)
+ if (ier /= 0) then
+ print*,'Error opening file: ',trim(CMTSOLUTION)
+ stop 'Error opening CMTSOLUTION file'
+ endif
! read source number isource
do isource=1,NSOURCES
- read(1,"(a)") string
+ ! initializes
+ yr = 0
+ da = 0
+ ho = -1
+ mi = -1
+ sec = -1.d0
+
+ ! gets header line
+ read(IIN,"(a256)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading header line in source ',isource
+ stop 'Error reading header line in station in CMTSOLUTION file'
+ endif
+
! skips empty lines
do while( len_trim(string) == 0 )
- read(1,"(a)") string
+ read(IIN,"(a256)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading header line in source ',isource
+ stop 'Error reading header line in station in CMTSOLUTION file'
+ endif
enddo
- ! read header with event information
- read(string,"(a4,i5,i3,i3,i3,i3,f6.2)") datasource,yr,mo,da,ho,mi,sec
+ ! debug
+ !print*,'line ----',string,'----'
+
+ ! reads header line with event information (assumes fixed format)
+ ! old line: read(string,"(a4,i5,i3,i3,i3,i3,f6.2)") datasource,yr,mo,da,ho,mi,sec
+
+ ! reads header line with event information (free format)
+ ! gets rid of the first datasource qualifyer string which can have variable length, like:
+ ! "PDE 2014 9 3 .."
+ ! " PDEQ2014 9 3 .."
+ ! " MLI 1971 1 1 .."
+ ! note: globalcmt.org solutions might have missing spaces after datasource qualifier
+ !
+ ! reads in year,month,day,hour,minutes,seconds
+ istart = 1
+ do itype = 1,6
+ ! determines where first number starts
+ do i = istart,len_trim(string)
+ if (is_numeric(string(i:i)) ) then
+ istart = i
+ exit
+ endif
+ enddo
+ if ( istart >= len_trim(string) ) stop 'Error determining datasource length in header line in CMTSOLUTION file'
+ if ( istart <= 1 ) stop 'Error determining datasource length in header line in CMTSOLUTION file'
+
+ ! determines end and length of number
+ iend = istart
+ do i = istart,len_trim(string)
+ if (itype /= 6 ) then
+ ! integer values
+ if (.not. is_numeric(string(i:i))) then
+ iend = i
+ exit
+ endif
+ else
+ ! seconds will have a digit number
+ ! digit numbers, e.g. 39.60, can contain '.'
+ if (.not. is_digit(string(i:i))) then
+ iend = i
+ exit
+ endif
+ endif
+ enddo
+ iend = iend-1
+ if ( iend >= len_trim(string) ) stop 'Error determining number length in header line in CMTSOLUTION file'
+ if ( iend < istart ) stop 'Error determining number with negative length in header line in CMTSOLUTION file'
+
+ ! debug
+ !print*,itype,'line ----',string(istart:iend),'----'
+
+ ! reads in event time information
+ select case( itype )
+ case( 1 )
+ ! year (as integer value)
+ read(string(istart:iend),*) yr
+ case( 2 )
+ ! month (as integer value)
+ read(string(istart:iend),*) mo
+ case( 3 )
+ ! day (as integer value)
+ read(string(istart:iend),*) da
+ case( 4 )
+ ! hour (as integer value)
+ read(string(istart:iend),*) ho
+ case( 5 )
+ ! minutes (as integer value)
+ read(string(istart:iend),*) mi
+ case( 6 )
+ ! seconds (as float value)
+ read(string(istart:iend),*) sec
+ end select
+
+ ! advances string
+ istart = iend + 1
+ enddo
+
+
+ ! checks time information
+ if (yr <= 0 .or. yr > 3000) then
+ print*, 'Error reading year: ',yr,' in source ',isource,'is invalid'
+ stop 'Error reading year out of header line in CMTSOLUTION file'
+ endif
+ if (mo < 1 .or. mo > 12) then
+ print*, 'Error reading month: ',mo,' in source ',isource,'is invalid'
+ stop 'Error reading month out of header line in CMTSOLUTION file'
+ endif
+ if (da < 1 .or. da > 31) then
+ print*, 'Error reading day: ',da,' in source ',isource,'is invalid'
+ stop 'Error reading day of header line in CMTSOLUTION file'
+ endif
+ if (ho < 0 .or. ho > 24) then
+ print*, 'Error reading hour: ',ho,' in source ',isource,'is invalid'
+ stop 'Error reading hour of header line in CMTSOLUTION file'
+ endif
+ if (mi < 0 .or. mi > 59) then
+ print*, 'Error reading minute: ',mi,' in source ',isource,'is invalid'
+ stop 'Error reading minute of header line in CMTSOLUTION file'
+ endif
+ if (sec < 0.0 .or. sec >= 60.0) then
+ print*, 'Error reading second: ',sec,' in source ',isource,'is invalid'
+ stop 'Error reading second of header line in CMTSOLUTION file'
+ endif
+
+ ! gets julian day number
jda=julian_day(yr,mo,da)
! ignore line with event name
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading event name in source ',isource
+ stop 'Error reading event name in station in CMTSOLUTION file'
+ endif
! read time shift
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading time shift in source ',isource
+ stop 'Error reading time shift in station in CMTSOLUTION file'
+ endif
+ !read(string(12:len_trim(string)),*) tshift_cmt(isource)
read(string(12:len_trim(string)),*) t_shift(isource)
! read half duration
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading half duration in source ',isource
+ stop 'Error reading half duration in station in CMTSOLUTION file'
+ endif
read(string(15:len_trim(string)),*) hdur(isource)
! read latitude
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading latitude in source ',isource
+ stop 'Error reading latitude in station in CMTSOLUTION file'
+ endif
read(string(10:len_trim(string)),*) lat(isource)
! read longitude
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading longitude in source ',isource
+ stop 'Error reading longitude in station in CMTSOLUTION file'
+ endif
read(string(11:len_trim(string)),*) long(isource)
! read depth
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading depth in source ',isource
+ stop 'Error reading depth in station in CMTSOLUTION file'
+ endif
read(string(7:len_trim(string)),*) depth(isource)
+ ! seismic moment tensor
+ ! CMTSOLUTION: components given in dyne-cm
! read Mrr
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading Mrr in source ',isource
+ stop 'Error reading Mrr in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(1,isource)
! read Mtt
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading Mtt in source ',isource
+ stop 'Error reading Mtt in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(2,isource)
! read Mpp
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading Mpp in source ',isource
+ stop 'Error reading Mpp in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(3,isource)
! read Mrt
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading Mrt in source ',isource
+ stop 'Error reading Mrt in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(4,isource)
! read Mrp
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading Mrp in source ',isource
+ stop 'Error reading Mrp in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(5,isource)
! read Mtp
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ print*, 'Error reading Mtp in source ',isource
+ stop 'Error reading Mtp in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(6,isource)
! checks half-duration
! null half-duration indicates a Heaviside
! replace with very short error function
- if( hdur(isource) < 5. * DT ) hdur(isource) = 5. * DT
+ if (hdur(isource) < 5. * DT ) hdur(isource) = 5. * DT
enddo
- close(1)
+ close(IIN)
! Sets tshift_cmt to zero to initiate the simulation!
if(NSOURCES == 1)then
@@ -156,8 +327,8 @@
endif
+ ! scales the moment tensor to Newton.m
!
- ! scale the moment tensor
! CMTSOLUTION file values are in dyne.cm
! 1 dyne is 1 gram * 1 cm / (1 second)^2
! 1 Newton is 1 kg * 1 m / (1 second)^2
@@ -166,8 +337,135 @@
!
moment_tensor(:,:) = moment_tensor(:,:) * 1.d-7
+ contains
+
+ !--------------------------------------------------------------
+
+ logical function is_numeric(char)
+
+ ! returns .true. if input character is a number
+
+ implicit none
+ character(len=1), intent(in) :: char
+
+ is_numeric = .false.
+
+ if ( index('0123456789', char) /= 0) then
+ is_numeric = .true.
+ endif
+
+ end function
+
+ !--------------------------------------------------------------
+
+ logical function is_digit(char)
+
+ ! returns .true. if input character is a number or a '.'
+
+ implicit none
+ character(len=1), intent(in) :: char
+
+ is_digit = .false.
+
+ if ( index('0123456789.', char) /= 0) then
+ is_digit = .true.
+ endif
+
+ end function
+
end subroutine get_cmt
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ double precision function get_cmt_scalar_moment(Mxx,Myy,Mzz,Mxy,Mxz,Myz)
+
+ ! calculates scalar moment (M0) in dyne-cm
+
+ implicit none
+
+ double precision, intent(in) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ ! local parameters
+ double precision :: scalar_moment,scaleM
+
+ ! scalar moment:
+ ! see equation (1.4) in P.G. Silver and T.H. Jordan, 1982,
+ ! "Optimal estiamtion of scalar seismic moment",
+ ! Geophys. J.R. astr. Soc., 70, 755 - 787
+ !
+ ! or see equation (5.91) in Dahlen & Tromp (1998)
+ !
+ ! moment tensor M is a symmetric 3x3 tensor, and has six independent components
+ !
+ ! the euclidean matrix norm is invariant under rotation.
+ ! thus, input can be:
+ ! Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ ! or
+ ! Mrr,Mtt,Mpp,Mrt,Mrp,Mtp
+ !
+ ! euclidean (or Frobenius) norm of a matrix: M0**2 = sum( Mij**2 )
+ scalar_moment = Mxx**2 + Myy**2 + Mzz**2 + 2.d0 * Mxy**2 + 2.d0 * Mxz**2 + 2.d0 * Myz**2
+
+ ! adds 1/2 to be coherent with double couple or point sources
+ scalar_moment = dsqrt(scalar_moment/2.0d0)
+
+ ! scale factor for the moment tensor
+ !
+ ! CMTSOLUTION file values are in Newton.m
+ ! 1 dyne is 1 gram * 1 cm / (1 second)^2
+ ! 1 Newton is 1 kg * 1 m / (1 second)^2
+ ! thus 1 Newton = 100,000 dynes
+ ! therefore 1 dyne.cm = 1e-7 Newton.m
+ ! and 1 Newton.m = 1e7 dyne.cm
+ scaleM = 1.d7
+
+ ! return value (in dyne-cm)
+ get_cmt_scalar_moment = scalar_moment * scaleM
+
+ end function
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ double precision function get_cmt_moment_magnitude(Mxx,Myy,Mzz,Mxy,Mxz,Myz)
+
+ ! calculates scalar moment (M0)
+
+ implicit none
+
+ double precision, intent(in) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ ! local parameters
+ double precision :: M0,Mw
+ double precision,external :: get_cmt_scalar_moment
+
+ ! scalar moment (in dyne-cm)
+ M0 = get_cmt_scalar_moment(Mxx,Myy,Mzz,Mxy,Mxz,Myz)
+
+ ! moment magnitude by Hanks & Kanamori, 1979
+ ! Mw = 2/3 log( M0 ) - 10.7 (dyne-cm)
+ !
+ ! alternative forms:
+ ! Mw = 2/3 ( log( M0 ) - 16.1 ) (N-m) "moment magnitude" by Hanks & Kanamori(1979) or "energy magnitude" by Kanamori (1977)
+ !
+ ! Aki & Richards ("Quantitative Seismology",2002):
+ ! Mw = 2/3 ( log( M0 ) - 9.1 ) (N-m)
+ !
+ ! conversion: dyne-cm = 10**-7 N-m
+ !
+ ! we follow here the USGS magnitude policy:
+ ! "All USGS statements of moment magnitude should use M = (log M0)/1.5-10.7
+ ! for converting from scalar moment M0 to moment magnitude. (..)"
+ ! see: http://earthquake.usgs.gov/aboutus/docs/020204mag_policy.php
+
+ Mw = 2.d0/3.d0 * log10( M0 ) - 10.7
+
+ ! return value
+ get_cmt_moment_magnitude = Mw
+
+ end function
+
! ------------------------------------------------------------------
integer function julian_day(yr,mo,da)
diff --git a/src/shared/get_force.f90 b/src/shared/get_force.f90
index b673138..0e9d24d 100644
--- a/src/shared/get_force.f90
+++ b/src/shared/get_force.f90
@@ -28,7 +28,7 @@
subroutine get_force(tshift_force,hdur,lat,long,depth,NSOURCES,min_tshift_force_original,factor_force_source, &
comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP)
- use constants
+ use constants,only: IIN,IN_DATA_FILES_PATH,MAX_STRING_LEN,TINYVAL,NUMBER_OF_SIMULTANEOUS_RUNS,mygroup
implicit none
@@ -42,13 +42,14 @@
double precision, dimension(NSOURCES), intent(out) :: comp_dir_vect_source_N
double precision, dimension(NSOURCES), intent(out) :: comp_dir_vect_source_Z_UP
-!--- local variables below
-
- integer isource,dummyval
- double precision t_shift(NSOURCES)
- double precision length
- character(len=7) dummy
- character(len=MAX_STRING_LEN) :: string, FORCESOLUTION,path_to_add
+ ! local variables below
+ integer :: isource,dummyval
+ double precision :: t_shift(NSOURCES)
+ double precision :: length
+ character(len=7) :: dummy
+ character(len=MAX_STRING_LEN) :: string
+ character(len=MAX_STRING_LEN) :: FORCESOLUTION,path_to_add
+ integer :: ier
! initializes
lat(:) = 0.d0
@@ -74,59 +75,63 @@
FORCESOLUTION = path_to_add(1:len_trim(path_to_add))//FORCESOLUTION(1:len_trim(FORCESOLUTION))
endif
- open(unit=1,file=trim(FORCESOLUTION),status='old',action='read')
+ open(unit=IIN,file=trim(FORCESOLUTION),status='old',action='read',iostat=ier)
+ if (ier /= 0) then
+ print*,'Error opening file: ',trim(FORCESOLUTION)
+ stop 'Error opening FORCESOLUTION file'
+ endif
! read source number isource
do isource=1,NSOURCES
- read(1,"(a)") string
+ read(IIN,"(a)") string
! skips empty lines
do while( len_trim(string) == 0 )
- read(1,"(a)") string
+ read(IIN,"(a)") string
enddo
! read header with event information
read(string,"(a6,i4)") dummy,dummyval
! read time shift
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(12:len_trim(string)),*) t_shift(isource)
! read half duration
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(6:len_trim(string)),*) hdur(isource)
! read latitude
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(10:len_trim(string)),*) lat(isource)
! read longitude
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(11:len_trim(string)),*) long(isource)
! read depth
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(7:len_trim(string)),*) depth(isource)
! read magnitude
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(21:len_trim(string)),*) factor_force_source(isource)
! read direction vector's East component
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(29:len_trim(string)),*) comp_dir_vect_source_E(isource)
! read direction vector's North component
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(29:len_trim(string)),*) comp_dir_vect_source_N(isource)
! read direction vector's vertical component
- read(1,"(a)") string
+ read(IIN,"(a)") string
read(string(32:len_trim(string)),*) comp_dir_vect_source_Z_UP(isource)
enddo
- close(1)
+ close(IIN)
! Sets tshift_force to zero to initiate the simulation!
if(NSOURCES == 1)then
diff --git a/src/shared/parallel.f90 b/src/shared/parallel.f90
index 2286a1b..81b287d 100644
--- a/src/shared/parallel.f90
+++ b/src/shared/parallel.f90
@@ -187,6 +187,25 @@ end module my_mpi
end subroutine bcast_all_dp
+
+!
+!----
+!
+
+ subroutine bcast_all_singledp(buffer)
+
+ use my_mpi
+
+ implicit none
+
+ double precision :: buffer
+
+ integer :: ier
+
+ call MPI_BCAST(buffer,1,MPI_DOUBLE_PRECISION,0,my_local_mpi_comm_world,ier)
+
+ end subroutine bcast_all_singledp
+
!
!----
!
diff --git a/src/shared/serial.f90 b/src/shared/serial.f90
index 2dc46ba..58b0fe4 100644
--- a/src/shared/serial.f90
+++ b/src/shared/serial.f90
@@ -126,6 +126,22 @@
!----
!
+ subroutine bcast_all_singledp(buffer)
+
+ use unused_mod
+
+ implicit none
+
+ double precision :: buffer
+
+ unused_dp = buffer
+
+ end subroutine bcast_all_singledp
+
+!
+!----
+!
+
subroutine bcast_all_r(buffer, count)
use unused_mod
diff --git a/src/specfem3D/locate_source.f90 b/src/specfem3D/locate_source.f90
index 23119ee..2fd62c3 100644
--- a/src/specfem3D/locate_source.f90
+++ b/src/specfem3D/locate_source.f90
@@ -35,24 +35,21 @@
DT,hdur,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
islice_selected_source,ispec_selected_source, &
xi_source,eta_source,gamma_source, &
- UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION,PRINT_SOURCE_TIME_FUNCTION, &
nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh, &
ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
use constants
- use specfem_par,only: USE_FORCE_POINT_SOURCE,USE_RICKER_TIME_FUNCTION,factor_force_source, &
- comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP
+ use specfem_par,only: USE_FORCE_POINT_SOURCE,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION, &
+ UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
+ factor_force_source,comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP
implicit none
- integer NPROC,UTM_PROJECTION_ZONE
+ integer NPROC
integer NSPEC_AB,NGLOB_AB,NSOURCES,NGNOD
- logical PRINT_SOURCE_TIME_FUNCTION
- logical SUPPRESS_UTM_PROJECTION
-
double precision DT
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -160,16 +157,41 @@
integer, dimension(NSOURCES) :: idomain
integer, dimension(NGATHER_SOURCES,0:NPROC-1) :: idomain_all
+ double precision, external :: get_cmt_scalar_moment
+ double precision, external :: get_cmt_moment_magnitude
+
!-----------------------------------------------------------------------------------
- ! read all the sources (note: each process reads the source file)
+ ! read all the sources
if (USE_FORCE_POINT_SOURCE) then
- call get_force(tshift_src,hdur,lat,long,depth,NSOURCES,min_tshift_src_original,factor_force_source, &
- comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP)
+ ! point forces
+ if (myrank == 0) then
+ ! only master process reads in FORCESOLUTION file
+ call get_force(tshift_src,hdur,lat,long,depth,NSOURCES,min_tshift_src_original,factor_force_source, &
+ comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP)
+ endif
+ ! broadcasts specific point force infos
+ call bcast_all_dp(factor_force_source,NSOURCES)
+ call bcast_all_dp(comp_dir_vect_source_E,NSOURCES)
+ call bcast_all_dp(comp_dir_vect_source_N,NSOURCES)
+ call bcast_all_dp(comp_dir_vect_source_Z_UP,NSOURCES)
else
- call get_cmt(yr,jda,ho,mi,sec,tshift_src,hdur,lat,long,depth,moment_tensor, &
- DT,NSOURCES,min_tshift_src_original)
+ ! CMT moment tensors
+ if (myrank == 0) then
+ ! only master process reads in CMTSOLUTION file
+ call get_cmt(yr,jda,ho,mi,sec,tshift_src,hdur,lat,long,depth,moment_tensor, &
+ DT,NSOURCES,min_tshift_src_original)
+ endif
+ ! broadcasts specific moment tensor infos
+ call bcast_all_dp(moment_tensor,6*NSOURCES)
endif
+ ! broadcasts general source information read on the master to the nodes
+ call bcast_all_dp(tshift_src,NSOURCES)
+ call bcast_all_dp(hdur,NSOURCES)
+ call bcast_all_dp(lat,NSOURCES)
+ call bcast_all_dp(long,NSOURCES)
+ call bcast_all_dp(depth,NSOURCES)
+ call bcast_all_singledp(min_tshift_src_original)
! define topology of the control element
call usual_hex_nodes(NGNOD,iaddx,iaddy,iaddz)
@@ -778,6 +800,13 @@
write(IMAIN,*)
endif
write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
+ write(IMAIN,*)
+ write(IMAIN,*) 'magnitude of the source:'
+ write(IMAIN,*) ' scalar moment M0 = ', &
+ get_cmt_scalar_moment(Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource)),' dyne-cm'
+ write(IMAIN,*) ' moment magnitude Mw = ', &
+ get_cmt_moment_magnitude(Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource))
+ write(IMAIN,*)
endif
write(IMAIN,*) ' time shift: ',tshift_src(isource),' seconds'
write(IMAIN,*)
diff --git a/src/specfem3D/prepare_timerun.F90 b/src/specfem3D/prepare_timerun.F90
index 270c610..03d1852 100644
--- a/src/specfem3D/prepare_timerun.F90
+++ b/src/specfem3D/prepare_timerun.F90
@@ -1194,6 +1194,9 @@
call flush_IMAIN()
endif
+ ! evaluates memory required
+ call memory_eval_gpu()
+
! prepares general fields on GPU
!§!§ JC JC here we will need to add GPU support for the new C-PML routines
call prepare_constants_device(Mesh_pointer, &
@@ -1364,6 +1367,160 @@
call flush_IMAIN()
endif
+ contains
+
+ subroutine memory_eval_gpu()
+
+ implicit none
+
+ ! local parameters
+ double precision :: memory_size
+ integer,parameter :: NGLL2 = 25
+ integer,parameter :: NGLL3 = 125
+ integer,parameter :: NGLL3_PADDED = 128
+
+ memory_size = 0.d0
+
+ ! add size of each set of arrays multiplied by the number of such arrays
+ ! d_hprime_xx,d_hprimewgll_xx
+ memory_size = memory_size + 2.d0 * NGLL2 * dble(CUSTOM_REAL)
+ ! padded xix,..gammaz
+ memory_size = memory_size + 9.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL)
+ ! padded kappav,muv
+ memory_size = memory_size + 2.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL)
+ ! ibool
+ memory_size = memory_size + NGLL3 * NSPEC_AB * dble(SIZE_INTEGER)
+ ! d_ibool_interfaces_ext_mesh
+ memory_size = memory_size + num_interfaces_ext_mesh * max_nibool_interfaces_ext_mesh * dble(SIZE_INTEGER)
+ ! ispec_is_inner
+ memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER)
+
+ if( STACEY_ABSORBING_CONDITIONS ) then
+ ! d_abs_boundary_ispec
+ memory_size = memory_size + num_abs_boundary_faces * dble(SIZE_INTEGER)
+ ! d_abs_boundary_ijk
+ memory_size = memory_size + num_abs_boundary_faces * 3.d0 * NGLL2 * dble(SIZE_INTEGER)
+ ! d_abs_boundary_normal
+ memory_size = memory_size + num_abs_boundary_faces * NDIM * NGLL2 * dble(CUSTOM_REAL)
+ ! d_abs_boundary_jacobian2Dw
+ memory_size = memory_size + num_abs_boundary_faces * NGLL2 * dble(CUSTOM_REAL)
+ endif
+
+ ! sources
+ ! d_sourcearrays
+ memory_size = memory_size + NGLL3 * NSOURCES * NDIM * dble(CUSTOM_REAL)
+ ! d_islice_selected_source,d_ispec_selected_source
+ memory_size = memory_size + 2.0 * NSOURCES * dble(SIZE_INTEGER)
+
+ ! receivers
+ !d_number_receiver_global
+ memory_size = memory_size + nrec_local * dble(SIZE_INTEGER)
+ ! d_ispec_selected_rec
+ memory_size = memory_size + nrec * dble(SIZE_INTEGER)
+
+ ! acoustic simulations
+ if( ACOUSTIC_SIMULATION ) then
+ ! d_potential_acoustic,d_potential_dot_acoustic,d_potential_dot_dot_acoustic
+ memory_size = memory_size + 3.d0 * NGLOB_AB * dble(CUSTOM_REAL)
+ ! d_rmass_acoustic
+ memory_size = memory_size + NGLOB_AB * dble(CUSTOM_REAL)
+ ! padded d_rhostore
+ memory_size = memory_size + NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL)
+ ! d_kappastore
+ memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL)
+ ! d_phase_ispec_inner_acoustic
+ memory_size = memory_size + 2.d0 * num_phase_ispec_acoustic * dble(SIZE_INTEGER)
+ ! d_ispec_is_acoustic
+ memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER)
+ endif
+
+ ! elastic simulations
+ if( ELASTIC_SIMULATION ) then
+ ! d_displ,..
+ memory_size = memory_size + 3.d0 * NDIM * NGLOB_AB * dble(CUSTOM_REAL)
+ ! d_send_accel_buffer
+ memory_size = memory_size + 3.d0 * num_interfaces_ext_mesh * max_nibool_interfaces_ext_mesh * dble(CUSTOM_REAL)
+ ! d_rmassx,..
+ memory_size = memory_size + 3.d0 * NGLOB_AB * dble(CUSTOM_REAL)
+ ! d_ispec_is_elastic
+ memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER)
+ ! d_phase_ispec_inner_elastic
+ memory_size = memory_size + 2.d0 * num_phase_ispec_elastic * dble(SIZE_INTEGER)
+ ! d_station_seismo_field
+ memory_size = memory_size + 3.d0 * NGLL3 * nrec_local * dble(CUSTOM_REAL)
+
+ if( STACEY_ABSORBING_CONDITIONS ) then
+ ! d_rho_vp,..
+ memory_size = memory_size + 2.d0 * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL)
+ endif
+ if( COMPUTE_AND_STORE_STRAIN ) then
+ ! d_epsilondev_xx,..
+ memory_size = memory_size + 5.d0 * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL)
+ endif
+ if( ATTENUATION ) then
+ ! d_R_xx,..
+ memory_size = memory_size + 5.d0 * size(R_xx) * dble(CUSTOM_REAL)
+ ! d_one_minus_sum_beta
+ memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL)
+ ! d_factor_common
+ memory_size = memory_size + N_SLS * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL)
+ ! alphaval,..
+ memory_size = memory_size + 3.d0 * N_SLS * dble(CUSTOM_REAL)
+ endif
+ if( ANISOTROPY ) then
+ ! padded d_c11store,..
+ memory_size = memory_size + 21.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL)
+ endif
+ if( APPROXIMATE_OCEAN_LOAD ) then
+ ! d_rmass_ocean_load
+ memory_size = memory_size + NGLOB_AB * dble(CUSTOM_REAL)
+ ! d_free_surface_normal
+ memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL)
+ endif
+ endif
+
+ ! noise simulations
+ if( NOISE_TOMOGRAPHY > 0 ) then
+ ! d_free_surface_ispec
+ memory_size = memory_size + num_free_surface_faces * dble(SIZE_INTEGER)
+ ! d_free_surface_ijk
+ memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(SIZE_INTEGER)
+ ! d_noise_surface_movie
+ memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL)
+ if( NOISE_TOMOGRAPHY == 1 ) then
+ ! d_noise_sourcearray
+ memory_size = memory_size + 3.d0 * NGLL3 * NSTEP * dble(CUSTOM_REAL)
+ endif
+ if( NOISE_TOMOGRAPHY > 1 ) then
+ ! d_normal_x_noise,..
+ memory_size = memory_size + 5.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL)
+ endif
+ if( NOISE_TOMOGRAPHY == 3 ) then
+ ! d_Sigma_kl
+ memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL)
+ endif
+ endif
+
+ if( GRAVITY ) then
+ ! d_minus_deriv_gravity,d_minus_g
+ memory_size = memory_size + 2.d0 * NGLOB_AB * dble(CUSTOM_REAL)
+ endif
+
+ ! poor estimate for kernel simulations...
+ if( SIMULATION_TYPE == 3 ) memory_size = 2.d0 * memory_size
+
+ ! user output
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) ' minimum memory requested : ', &
+ memory_size / 1024. / 1024., &
+ 'MB per process'
+ write(IMAIN,*)
+ call flush_IMAIN()
+ endif
+
+ end subroutine memory_eval_gpu
+
end subroutine prepare_timerun_GPU
!
diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90
index 44e0184..50d9655 100644
--- a/src/specfem3D/setup_sources_receivers.f90
+++ b/src/specfem3D/setup_sources_receivers.f90
@@ -126,7 +126,6 @@
DT,hdur,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
islice_selected_source,ispec_selected_source, &
xi_source,eta_source,gamma_source, &
- UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION,PRINT_SOURCE_TIME_FUNCTION, &
nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh,&
ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
More information about the CIG-COMMITS
mailing list