[cig-commits] r12654 - seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Aug 16 09:22:55 PDT 2008
Author: dkomati1
Date: 2008-08-16 09:22:54 -0700 (Sat, 16 Aug 2008)
New Revision: 12654
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
Log:
modified the implementation of attenuation to improve performance: inlined and simplified the calls to Brian Savage's routines
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90 2008-08-16 15:29:47 UTC (rev 12653)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90 2008-08-16 16:22:54 UTC (rev 12654)
@@ -887,6 +887,8 @@
end subroutine set_attenuation_regions_1D
+!------------------------------------------------------------------------
+
subroutine get_attenuation_index(iflag, radius, index, inner_core, AM_V)
implicit none
@@ -972,7 +974,7 @@
endif
-! Clamp regions
+! clamp regions
if(index < AM_V%Qrmin(iregion)) index = AM_V%Qrmin(iregion)
if(index > AM_V%Qrmax(iregion)) index = AM_V%Qrmax(iregion)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90 2008-08-16 15:29:47 UTC (rev 12653)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90 2008-08-16 16:22:54 UTC (rev 12654)
@@ -83,7 +83,7 @@
real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
- integer iregion_selected
+ integer iregion_selected,iregion
! for attenuation
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
@@ -256,9 +256,24 @@
! precompute terms for attenuation if needed
if(ATTENUATION_VAL) then
- radius_cr = xstore(ibool(i,j,k,ispec))
- call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
+!! DK DK do not call Brian's routine anymore, to improve performance I rewrote it differently
+ radius_cr = xstore(ibool(i,j,k,ispec))
+!!! call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
+ iregion_selected = nint(radius_cr * TABLE_ATTENUATION)
+ if(idoubling(ispec) == IFLAG_MANTLE_NORMAL) then
+ iregion = IREGION_ATTENUATION_CMB_670
+ else if(idoubling(ispec) == IFLAG_670_220) then
+ iregion = IREGION_ATTENUATION_670_220
+ else if(idoubling(ispec) == IFLAG_220_80) then
+ iregion = IREGION_ATTENUATION_220_80
+ else if(idoubling(ispec) == IFLAG_CRUST .or. idoubling(ispec) == IFLAG_80_MOHO) then
+ iregion = IREGION_ATTENUATION_80_SURFACE
+ endif
+! clamp regions
+ if(iregion_selected < AM_V%Qrmin(iregion)) iregion_selected = AM_V%Qrmin(iregion)
+ if(iregion_selected > AM_V%Qrmax(iregion)) iregion_selected = AM_V%Qrmax(iregion)
+!! DK DK do not call Brian's routine anymore, to improve performance I rewrote it differently
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
minus_sum_beta = one_minus_sum_beta_use - 1.0
endif
@@ -594,35 +609,35 @@
enddo
enddo
-! update memory variables based upon the Runge-Kutta scheme
-! convention for attenuation
+! update memory variables based upon a Runge-Kutta scheme.
+! convention for attenuation:
! term in xx = 1
! term in yy = 2
! term in xy = 3
! term in xz = 4
! term in yz = 5
! term in zz not computed since zero trace
-
- if(ATTENUATION_VAL) then
-
-! use Runge-Kutta scheme to march in time
- do i_sls = 1,N_SLS
- do i_memory = 1,5
-
+ if(ATTENUATION_VAL) then
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ do i_sls = 1,N_SLS
+ do i_memory = 1,5
! get coefficients for that standard linear solid
! IMPROVE we use mu_v here even if there is some anisotropy
! IMPROVE we should probably use an average value instead
+ R_memory(i_memory,i_sls,i,j,k,ispec) = alphaval(i_sls) * &
+ R_memory(i_memory,i_sls,i,j,k,ispec) + &
+ factor_common(i_sls,1,1,1,iregion_selected) * muvstore(i,j,k,ispec) * &
+ (betaval(i_sls) * epsilondev(i_memory,i,j,k,ispec) + &
+ gammaval(i_sls) * epsilondev_loc(i_memory,i,j,k))
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
- R_memory(i_memory,i_sls,:,:,:,ispec) = alphaval(i_sls) * &
- R_memory(i_memory,i_sls,:,:,:,ispec) + &
- factor_common(i_sls,1,1,1,iregion_selected) * muvstore(:,:,:,ispec) * &
- (betaval(i_sls) * epsilondev(i_memory,:,:,:,ispec) + &
- gammaval(i_sls) * epsilondev_loc(i_memory,:,:,:))
- enddo
- enddo
-
- endif
-
! save deviatoric strain for Runge-Kutta scheme
if(COMPUTE_AND_STORE_STRAIN) epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90 2008-08-16 15:29:47 UTC (rev 12653)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90 2008-08-16 16:22:54 UTC (rev 12654)
@@ -25,7 +25,7 @@
!
!=====================================================================
- subroutine compute_forces_inner_core(displ,accel,xstore, &
+ subroutine compute_forces_inner_core(displ,accel, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
@@ -33,7 +33,7 @@
kappavstore,muvstore,ibool,idoubling, &
R_memory,epsilondev,&
one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+ vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN)
implicit none
@@ -43,28 +43,9 @@
! done for performance only using static allocation to allow for loop unrolling
include "values_from_mesher.h"
-! attenuation_model_variables
- type attenuation_model_variables
- sequence
- double precision min_period, max_period
- double precision :: QT_c_source ! Source Frequency
- double precision, dimension(N_SLS) :: Qtau_s ! tau_sigma
- double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
- double precision, dimension(:), pointer :: Qr ! Radius
- integer, dimension(:), pointer :: interval_Q ! Steps
- double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
- double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
- double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
- double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
- double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
- integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
- integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
- integer :: Qn ! Number of points
- end type attenuation_model_variables
+! same attenuation everywhere in the inner core therefore no need to use Brian's routines
+ integer, parameter :: iregion_selected = 1
- type (attenuation_model_variables) AM_V
-! attenuation_model_variables
-
! for forward or backward simulations
logical COMPUTE_AND_STORE_STRAIN
@@ -133,12 +114,6 @@
real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
-! for gravity
- real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: xstore
- integer iregion_selected
-
- real(kind=CUSTOM_REAL) radius_cr
-
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -238,8 +213,9 @@
endif
if(ATTENUATION_VAL) then
- radius_cr = xstore(ibool(i,j,k,ispec))
- call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
+! same attenuation everywhere in the inner core therefore no need to use Brian's routines
+!!!!! radius_cr = xstore(ibool(i,j,k,ispec))
+!!!!! call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
minus_sum_beta = one_minus_sum_beta(1,1,1,iregion_selected) - 1.0
endif ! ATTENUATION_VAL
@@ -362,34 +338,35 @@
enddo
enddo
-! use Runge-Kutta scheme to march memory variables in time
-! convention for attenuation
+! update memory variables based upon a Runge-Kutta scheme.
+! convention for attenuation:
! term in xx = 1
! term in yy = 2
! term in xy = 3
! term in xz = 4
! term in yz = 5
! term in zz not computed since zero trace
-
if(ATTENUATION_VAL) then
-
- do i_sls = 1,N_SLS
- do i_memory = 1,5
- R_memory(i_memory,i_sls,:,:,:,ispec) = &
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ do i_sls = 1,N_SLS
+ do i_memory = 1,5
+ R_memory(i_memory,i_sls,i,j,k,ispec) = &
alphaval(i_sls) * &
- R_memory(i_memory,i_sls,:,:,:,ispec) + muvstore(:,:,:,ispec) * &
+ R_memory(i_memory,i_sls,i,j,k,ispec) + muvstore(i,j,k,ispec) * &
factor_common(i_sls,1,1,1,iregion_selected) * &
(betaval(i_sls) * &
- epsilondev(i_memory,:,:,:,ispec) + gammaval(i_sls) * epsilondev_loc(i_memory,:,:,:))
+ epsilondev(i_memory,i,j,k,ispec) + gammaval(i_sls) * epsilondev_loc(i_memory,i,j,k))
+ enddo
+ enddo
enddo
- enddo
+ enddo
+ enddo
+ endif
- endif
-
- if (COMPUTE_AND_STORE_STRAIN) then
! save deviatoric strain for Runge-Kutta scheme
- epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
- endif
+ if(COMPUTE_AND_STORE_STRAIN) epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
endif ! end test to exclude fictitious elements in central cube
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-08-16 15:29:47 UTC (rev 12653)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-08-16 16:22:54 UTC (rev 12654)
@@ -1452,7 +1452,6 @@
! W. H. Freeman, (1980), second edition, sections 5.5 and 5.5.2, eq. (5.81) p. 170
! rescale in crust and mantle
-
do ispec = 1,NSPEC_CRUST_MANTLE
do k=1,NGLLZ
do j=1,NGLLY
@@ -1503,7 +1502,6 @@
enddo ! END DO CRUST MANTLE
! rescale in inner core
-
do ispec = 1,NSPEC_INNER_CORE
do k=1,NGLLZ
do j=1,NGLLY
@@ -2157,7 +2155,7 @@
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
- call compute_forces_inner_core(displ_inner_core,accel_inner_core,xstore_inner_core, &
+ call compute_forces_inner_core(displ_inner_core,accel_inner_core, &
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
@@ -2168,7 +2166,7 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN)
#ifdef USE_MPI
More information about the cig-commits
mailing list