[cig-commits] r22584 - in seismo/3D/SPECFEM3D_GLOBE/trunk: setup src/specfem3D utils
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Jul 13 07:08:06 PDT 2013
Author: dkomati1
Date: 2013-07-13 07:08:06 -0700 (Sat, 13 Jul 2013)
New Revision: 22584
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/expand_compute_strain_product.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_kernels.f90
Log:
started to vectorize the code using the FORCE_VECTORIZATION flag;
also moved rotate_kernels_dble() from compute_kernels.F90 to save_kernels.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2013-07-13 14:03:38 UTC (rev 22583)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2013-07-13 14:08:06 UTC (rev 22584)
@@ -349,6 +349,7 @@
! Deville routines optimized for NGLLX = NGLLY = NGLLZ = 5
integer, parameter :: m1 = NGLLX, m2 = NGLLX * NGLLY
+ integer, parameter :: NGLLCUBE = NGLLX * NGLLY * NGLLZ
! mid-points inside a GLL element
integer, parameter :: MIDX = (NGLLX+1)/2
@@ -360,6 +361,7 @@
! a few useful constants
double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
+ double precision, parameter :: ONE_HALF = HALF, ONE_FOURTH = 0.25d0, ONE_EIGHTH = 0.125d0, ONE_SIXTEENTH = 0.0625d0
real(kind=CUSTOM_REAL), parameter :: &
ONE_THIRD = 1._CUSTOM_REAL/3._CUSTOM_REAL, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-07-13 14:03:38 UTC (rev 22583)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-07-13 14:08:06 UTC (rev 22584)
@@ -237,6 +237,9 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#endif
! ****************************************************
! big loop over all spectral elements in the solid
@@ -287,6 +290,14 @@
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ iglob1 = ibool(ijk,1,1,ispec)
+ dummyx_loc(ijk,1,1) = displ_crust_mantle(1,iglob1)
+ dummyy_loc(ijk,1,1) = displ_crust_mantle(2,iglob1)
+ dummyz_loc(ijk,1,1) = displ_crust_mantle(3,iglob1)
+ enddo
+#else
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -297,6 +308,7 @@
enddo
enddo
enddo
+#endif
do j=1,m2
do i=1,m1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-07-13 14:03:38 UTC (rev 22583)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-07-13 14:08:06 UTC (rev 22584)
@@ -252,6 +252,10 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#endif
+
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -300,7 +304,14 @@
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
-
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ iglob1 = ibool(ijk,1,1,ispec)
+ dummyx_loc(ijk,1,1) = displ_inner_core(1,iglob1)
+ dummyy_loc(ijk,1,1) = displ_inner_core(2,iglob1)
+ dummyz_loc(ijk,1,1) = displ_inner_core(3,iglob1)
+ enddo
+#else
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -311,6 +322,7 @@
enddo
enddo
enddo
+#endif
do j=1,m2
do i=1,m1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90 2013-07-13 14:03:38 UTC (rev 22583)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90 2013-07-13 14:08:06 UTC (rev 22584)
@@ -134,15 +134,15 @@
beta_kl_crust_mantle(i,j,k,ispec) = beta_kl_crust_mantle(i,j,k,ispec) &
+ deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
+ (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
- + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
+ + 2.d0 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
epsilondev_loc(5)*b_epsilondev_loc(5)) )
! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
! note: multiplication with kappa(x) will be done after the time loop
alpha_kl_crust_mantle(i,j,k,ispec) = alpha_kl_crust_mantle(i,j,k,ispec) &
- + deltat * (9 * eps_trace_over_3_loc_matrix(i,j,k) &
- * b_eps_trace_over_3_loc_matrix(i,j,k))
+ + deltat * (9.d0 * eps_trace_over_3_loc_matrix(i,j,k) &
+ * b_eps_trace_over_3_loc_matrix(i,j,k))
endif
@@ -311,7 +311,6 @@
end subroutine compute_kernels_outer_core
-
!
!-------------------------------------------------------------------------------------------------
!
@@ -350,14 +349,18 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: eps_trace_over_3_inner_core
! local parameters
- real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
- real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_matrix,b_epsilondev_loc_matrix
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_matrix,&
b_eps_trace_over_3_loc_matrix
- integer :: i,j,k,ispec,iglob
+ integer :: ispec,iglob
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#else
+ integer :: i,j,k
+#endif
+
! inner_core
do ispec = 1, NSPEC_INNER_CORE
@@ -380,6 +383,28 @@
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,&
b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
+#ifdef FORCE_VECTORIZATION
+ do ijk = 1,NGLLCUBE
+ iglob = ibool_inner_core(ijk,1,1,ispec)
+
+ rho_kl_inner_core(ijk,1,1,ispec) = rho_kl_inner_core(ijk,1,1,ispec) &
+ + deltat * (accel_inner_core(1,iglob) * b_displ_inner_core(1,iglob) &
+ + accel_inner_core(2,iglob) * b_displ_inner_core(2,iglob) &
+ + accel_inner_core(3,iglob) * b_displ_inner_core(3,iglob))
+
+ beta_kl_inner_core(ijk,1,1,ispec) = beta_kl_inner_core(ijk,1,1,ispec) &
+ + deltat * (epsilondev_loc_matrix(1,ijk,1,1)*b_epsilondev_loc_matrix(1,ijk,1,1) &
+ + epsilondev_loc_matrix(2,ijk,1,1)*b_epsilondev_loc_matrix(2,ijk,1,1) &
+ + (epsilondev_loc_matrix(1,ijk,1,1)+epsilondev_loc_matrix(2,ijk,1,1)) * &
+ (b_epsilondev_loc_matrix(1,ijk,1,1)+b_epsilondev_loc_matrix(2,ijk,1,1)) &
+ + 2.d0 * (epsilondev_loc_matrix(3,ijk,1,1)*b_epsilondev_loc_matrix(3,ijk,1,1) &
+ + epsilondev_loc_matrix(4,ijk,1,1)*b_epsilondev_loc_matrix(4,ijk,1,1) &
+ + epsilondev_loc_matrix(5,ijk,1,1)*b_epsilondev_loc_matrix(5,ijk,1,1)))
+
+ alpha_kl_inner_core(ijk,1,1,ispec) = alpha_kl_inner_core(ijk,1,1,ispec) &
+ + deltat * (9.d0 * eps_trace_over_3_loc_matrix(ijk,1,1) * b_eps_trace_over_3_loc_matrix(ijk,1,1))
+ enddo
+#else
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
@@ -388,22 +413,24 @@
rho_kl_inner_core(i,j,k,ispec) = rho_kl_inner_core(i,j,k,ispec) &
+ deltat * (accel_inner_core(1,iglob) * b_displ_inner_core(1,iglob) &
+ accel_inner_core(2,iglob) * b_displ_inner_core(2,iglob) &
- + accel_inner_core(3,iglob) * b_displ_inner_core(3,iglob) )
+ + accel_inner_core(3,iglob) * b_displ_inner_core(3,iglob))
- epsilondev_loc(:) = epsilondev_loc_matrix(:,i,j,k)
- b_epsilondev_loc(:) = b_epsilondev_loc_matrix(:,i,j,k)
-
beta_kl_inner_core(i,j,k,ispec) = beta_kl_inner_core(i,j,k,ispec) &
- + deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
- + (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
- + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) &
- + epsilondev_loc(5)*b_epsilondev_loc(5)) )
+ + deltat * (epsilondev_loc_matrix(1,i,j,k)*b_epsilondev_loc_matrix(1,i,j,k) &
+ + epsilondev_loc_matrix(2,i,j,k)*b_epsilondev_loc_matrix(2,i,j,k) &
+ + (epsilondev_loc_matrix(1,i,j,k)+epsilondev_loc_matrix(2,i,j,k)) * &
+ (b_epsilondev_loc_matrix(1,i,j,k)+b_epsilondev_loc_matrix(2,i,j,k)) &
+ + 2.d0 * (epsilondev_loc_matrix(3,i,j,k)*b_epsilondev_loc_matrix(3,i,j,k) &
+ + epsilondev_loc_matrix(4,i,j,k)*b_epsilondev_loc_matrix(4,i,j,k) &
+ + epsilondev_loc_matrix(5,i,j,k)*b_epsilondev_loc_matrix(5,i,j,k)))
alpha_kl_inner_core(i,j,k,ispec) = alpha_kl_inner_core(i,j,k,ispec) &
- + deltat * (9 * eps_trace_over_3_loc_matrix(i,j,k) * b_eps_trace_over_3_loc_matrix(i,j,k))
+ + deltat * (9.d0 * eps_trace_over_3_loc_matrix(i,j,k) * b_eps_trace_over_3_loc_matrix(i,j,k))
enddo
enddo
enddo
+#endif
+
enddo
end subroutine compute_kernels_inner_core
@@ -424,7 +451,6 @@
! prod(7)=eps22*eps22 -> c22, prod(8)=eps22eps33 -> c23, prod(9)=eps22eps23 -> c24, ...
! prod(19)=eps13*eps13 -> c55, prod(20)=eps13eps12 -> c56, prod(21)=eps12eps12 -> c66
! This then gives how the 21 kernels are organized
- ! For crust_mantle
implicit none
include "constants.h"
@@ -438,11 +464,11 @@
! Building of the local matrix of the strain tensor
! for the adjoint field and the regular backward field
- eps(1:2)=epsdev(1:2)+eps_trace_over_3 !eps11 and eps22
- eps(3)=-(eps(1)+eps(2))+3*eps_trace_over_3 !eps33
- eps(4)=epsdev(5) !eps23
- eps(5)=epsdev(4) !eps13
- eps(6)=epsdev(3) !eps12
+ eps(1:2)=epsdev(1:2)+eps_trace_over_3 ! eps11 and eps22
+ eps(3)=-(eps(1)+eps(2))+3*eps_trace_over_3 ! eps33
+ eps(4)=epsdev(5) ! eps23
+ eps(5)=epsdev(4) ! eps13
+ eps(6)=epsdev(3) ! eps12
b_eps(1:2)=b_epsdev(1:2)+b_eps_trace_over_3
b_eps(3)=-(b_eps(1)+b_eps(2))+3*b_eps_trace_over_3
@@ -470,383 +496,6 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine rotate_kernels_dble(cij_kl,cij_kll,theta_in,phi_in)
-
-! Purpose : compute the kernels in r,theta,phi (cij_kll)
-! from the kernels in x,y,z (cij_kl) (x,y,z <-> r,theta,phi)
-! At r,theta,phi fixed
-! theta and phi are in radians
-
-! Coeff from Min's routine rotate_anisotropic_tensor
-! with the help of Collect[Expand[cij],{dij}] in Mathematica
-
-! Definition of the output array cij_kll :
-! cij_kll(1) = C11 ; cij_kll(2) = C12 ; cij_kll(3) = C13
-! cij_kll(4) = C14 ; cij_kll(5) = C15 ; cij_kll(6) = C16
-! cij_kll(7) = C22 ; cij_kll(8) = C23 ; cij_kll(9) = C24
-! cij_kll(10) = C25 ; cij_kll(11) = C26 ; cij_kll(12) = C33
-! cij_kll(13) = C34 ; cij_kll(14) = C35 ; cij_kll(15) = C36
-! cij_kll(16) = C44 ; cij_kll(17) = C45 ; cij_kll(18) = C46
-! cij_kll(19) = C55 ; cij_kll(20) = C56 ; cij_kll(21) = C66
-! where the Cij (Voigt's notation) are defined as function of
-! the components of the elastic tensor in spherical coordinates
-! by eq. (A.1) of Chen & Tromp, GJI 168 (2007)
-
- implicit none
- include "constants.h"
-
- real(kind=CUSTOM_REAL), intent(in) :: theta_in,phi_in
-
- real(kind=CUSTOM_REAL), dimension(21), intent(in) :: cij_kl
- real(kind=CUSTOM_REAL), dimension(21), intent(out) :: cij_kll
-
- double precision :: theta,phi
- double precision :: costheta,sintheta,cosphi,sinphi
- double precision :: costhetasq,sinthetasq,cosphisq,sinphisq
- double precision :: costwotheta,sintwotheta,costwophi,sintwophi
- double precision :: cosfourtheta,sinfourtheta,cosfourphi,sinfourphi
- double precision :: costhetafour,sinthetafour,cosphifour,sinphifour
- double precision :: sintwophisq,sintwothetasq
- double precision :: costhreetheta,sinthreetheta,costhreephi,sinthreephi
-
- if (CUSTOM_REAL == SIZE_REAL) then
- theta = dble(theta_in)
- phi = dble(phi_in)
- else
- theta = theta_in
- phi = phi_in
- endif
-
- costheta = dcos(theta)
- sintheta = dsin(theta)
- cosphi = dcos(phi)
- sinphi = dsin(phi)
-
- costhetasq = costheta * costheta
- sinthetasq = sintheta * sintheta
- cosphisq = cosphi * cosphi
- sinphisq = sinphi * sinphi
-
- costhetafour = costhetasq * costhetasq
- sinthetafour = sinthetasq * sinthetasq
- cosphifour = cosphisq * cosphisq
- sinphifour = sinphisq * sinphisq
-
- costwotheta = dcos(2.d0*theta)
- sintwotheta = dsin(2.d0*theta)
- costwophi = dcos(2.d0*phi)
- sintwophi = dsin(2.d0*phi)
-
- costhreetheta=dcos(3.d0*theta)
- sinthreetheta=dsin(3.d0*theta)
- costhreephi=dcos(3.d0*phi)
- sinthreephi=dsin(3.d0*phi)
-
- cosfourtheta = dcos(4.d0*theta)
- sinfourtheta = dsin(4.d0*theta)
- cosfourphi = dcos(4.d0*phi)
- sinfourphi = dsin(4.d0*phi)
- sintwothetasq = sintwotheta * sintwotheta
- sintwophisq = sintwophi * sintwophi
-
-
- cij_kll(1) = 1.d0/16.d0* (cij_kl(16) - cij_kl(16)* costwophi + &
- 16.d0* cosphi*cosphisq* costhetafour* (cij_kl(1)* cosphi + cij_kl(6)* sinphi) + &
- 2.d0* (cij_kl(15) + cij_kl(17))* sintwophi* sintwothetasq - &
- 2.d0* (cij_kl(16)* cosfourtheta* sinphisq + &
- 2.d0* costhetafour* (-4* cij_kl(7)* sinphifour - &
- (cij_kl(2) + cij_kl(21))* sintwophisq) + &
- 8.d0* cij_kl(5)* cosphi*cosphisq* costheta*costhetasq* sintheta - &
- 8.d0* cij_kl(8)* costhetasq* sinphisq* sinthetasq - &
- 8.d0* cij_kl(12)* sinthetafour + &
- 8.d0* cosphisq* costhetasq* sintheta* ((cij_kl(4) + &
- cij_kl(20))* costheta* sinphi - &
- (cij_kl(3) + cij_kl(19))*sintheta) + &
- 8.d0* cosphi* costheta* (-cij_kl(11)* costheta*costhetasq* &
- sinphi*sinphisq + (cij_kl(10) + cij_kl(18))* costhetasq* sinphisq* sintheta + &
- cij_kl(14)* sintheta*sinthetasq) + 2.d0* sinphi* (cij_kl(13) + &
- cij_kl(9)* sinphisq)* sintwotheta + &
- sinphi* (-cij_kl(13) + cij_kl(9)* sinphisq)* sinfourtheta))
-
- cij_kll(2) = 1.d0/4.d0* (costhetasq* (cij_kl(1) + 3.d0* cij_kl(2) + cij_kl(7) - &
- cij_kl(21) + (-cij_kl(1) + cij_kl(2) - cij_kl(7) + &
- cij_kl(21))* cosfourphi + (-cij_kl(6) + cij_kl(11))* sinfourphi) + &
- 4.d0* (cij_kl(8)* cosphisq - cij_kl(15)* cosphi* sinphi + &
- cij_kl(3)* sinphisq)* sinthetasq - &
- 2.d0* (cij_kl(10)* cosphisq*cosphi + &
- (cij_kl(9) - cij_kl(20))* cosphisq* sinphi + &
- (cij_kl(5) - cij_kl(18))* cosphi* sinphisq + &
- cij_kl(4)* sinphisq*sinphi)* sintwotheta)
-
- cij_kll(3) = 1.d0/8.d0* (sintwophi* (3.d0* cij_kl(15) - cij_kl(17) + &
- 4.d0* (cij_kl(2) + cij_kl(21))* costhetasq* sintwophi* sinthetasq) + &
- 4.d0* cij_kl(12)* sintwothetasq + 4.d0* cij_kl(1)* cosphifour* sintwothetasq + &
- 2.d0* cosphi*cosphisq* (8.d0* cij_kl(6)* costhetasq* sinphi* sinthetasq + &
- cij_kl(5)* sinfourtheta) + 2.d0* cosphisq* (3.d0* cij_kl(3) - cij_kl(19) + &
- (cij_kl(3) + cij_kl(19))* cosfourtheta + &
- (cij_kl(4) + cij_kl(20))* sinphi* sinfourtheta) + &
- 2.d0* sinphi* (sinphi* (3.d0* cij_kl(8) - &
- cij_kl(16) + (cij_kl(8) + cij_kl(16))* cosfourtheta + &
- 2.d0* cij_kl(7)* sinphisq* sintwothetasq)+ &
- (-cij_kl(13) + cij_kl(9)* sinphisq)* sinfourtheta)+ &
- 2.d0* cosphi* ((cij_kl(15) + cij_kl(17))* cosfourtheta* sinphi + &
- 8.d0* cij_kl(11)* costhetasq* sinphi*sinphisq* sinthetasq + &
- (-cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq)*sinfourtheta))
-
- cij_kll(4) = 1.d0/8.d0* (cosphi* costheta *(5.d0* cij_kl(4) - &
- cij_kl(9) + 4.d0* cij_kl(13) - &
- 3.d0* cij_kl(20) + (cij_kl(4) + 3.d0* cij_kl(9) - &
- 4.d0* cij_kl(13) + cij_kl(20))* costwotheta) + &
- 1.d0/2.d0* (cij_kl(4) - cij_kl(9) + &
- cij_kl(20))* costhreephi * (costheta + 3.d0* costhreetheta) - &
- costheta* (-cij_kl(5) + 5.d0* cij_kl(10) + &
- 4.d0* cij_kl(14) - 3.d0* cij_kl(18) + &
- (3.d0* cij_kl(5) + cij_kl(10) - &
- 4.d0* cij_kl(14) + cij_kl(18))* costwotheta)* sinphi - &
- 1.d0/2.d0* (cij_kl(5) - cij_kl(10) - cij_kl(18))* (costheta + &
- 3.d0* costhreetheta)* sinthreephi + &
- 4.d0* (cij_kl(6) - cij_kl(11))* cosfourphi* costhetasq* sintheta - &
- 4.d0* (cij_kl(1) + cij_kl(3) - cij_kl(7) - cij_kl(8) + cij_kl(16) - cij_kl(19) + &
- (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + &
- cij_kl(16) - cij_kl(19))* costwotheta)* sintwophi* sintheta - &
- 4.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
- cij_kl(21))* costhetasq* sinfourphi* sintheta + &
- costwophi* ((cij_kl(6) + cij_kl(11) + 6.d0* cij_kl(15) - &
- 2.d0* cij_kl(17))* sintheta + &
- (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + cij_kl(17)))* sinthreetheta))
-
- cij_kll(5) = 1.d0/4.d0* (2.d0* (cij_kl(4) + &
- cij_kl(20))* cosphisq* (costwotheta + cosfourtheta)* sinphi + &
- 2.d0* cij_kl(9)* (costwotheta + cosfourtheta)* sinphi*sinphisq + &
- 16.d0* cij_kl(1)* cosphifour* costheta*costhetasq* sintheta + &
- 4.d0* costheta*costhetasq* (-2.d0* cij_kl(8)* sinphisq + &
- 4.d0* cij_kl(7)* sinphifour + &
- (cij_kl(2) + cij_kl(21))* sintwophisq)* sintheta + &
- 4.d0* cij_kl(13)* (1.d0 + 2.d0* costwotheta)* sinphi* sinthetasq + &
- 8.d0* costheta* (-2.d0* cij_kl(12) + cij_kl(8)* sinphisq)* sintheta*sinthetasq + &
- 2.d0* cosphi*cosphisq* (cij_kl(5)* (costwotheta + cosfourtheta) + &
- 8.d0* cij_kl(6)* costheta*costhetasq* sinphi* sintheta) + &
- 2.d0* cosphi* (cosfourtheta* (-cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq) + &
- costwotheta* (cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq) + &
- 8.d0* cij_kl(11)* costheta*costhetasq* sinphi*sinphisq* sintheta) - &
- (cij_kl(3) + cij_kl(16) + cij_kl(19) + &
- (cij_kl(3) - cij_kl(16) + cij_kl(19))* costwophi + &
- (cij_kl(15) + cij_kl(17))* sintwophi)* sinfourtheta)
-
- cij_kll(6) = 1.d0/2.d0* costheta*costhetasq* ((cij_kl(6) + cij_kl(11))* costwophi + &
- (cij_kl(6) - cij_kl(11))* cosfourphi + 2.d0* (-cij_kl(1) + cij_kl(7))* sintwophi + &
- (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* sinfourphi) + &
- 1.d0/4.d0* costhetasq* (-(cij_kl(4) + 3* cij_kl(9) + cij_kl(20))* cosphi - &
- 3.d0* (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi + &
- (3.d0* cij_kl(5) + cij_kl(10) + cij_kl(18))* sinphi + &
- 3.d0* (cij_kl(5) - cij_kl(10) - cij_kl(18))* sinthreephi)* sintheta + &
- costheta* ((cij_kl(15) + cij_kl(17))* costwophi + &
- (-cij_kl(3) + cij_kl(8) + cij_kl(16) - cij_kl(19))* sintwophi)* sinthetasq + &
- (-cij_kl(13)* cosphi + cij_kl(14)* sinphi)* sintheta*sinthetasq
-
- cij_kll(7) = cij_kl(7)* cosphifour - cij_kl(11)* cosphi*cosphisq* sinphi + &
- (cij_kl(2) + cij_kl(21))* cosphisq* sinphisq - &
- cij_kl(6)* cosphi* sinphi*sinphisq + &
- cij_kl(1)* sinphifour
-
- cij_kll(8) = 1.d0/2.d0* (2.d0* costhetasq* sinphi* (-cij_kl(15)* cosphi + &
- cij_kl(3)* sinphi) + 2.d0* cij_kl(2)* cosphifour* sinthetasq + &
- (2.d0* cij_kl(2)* sinphifour + &
- (cij_kl(1) + cij_kl(7) - cij_kl(21))* sintwophisq)* sinthetasq + &
- cij_kl(4)* sinphi*sinphisq* sintwotheta + &
- cosphi*cosphisq* (2.d0* (-cij_kl(6) + cij_kl(11))* sinphi* sinthetasq + &
- cij_kl(10)* sintwotheta) + cosphi* sinphisq* (2.d0* (cij_kl(6) - &
- cij_kl(11))* sinphi* sinthetasq + &
- (cij_kl(5) - cij_kl(18))* sintwotheta) + &
- cosphisq* (2.d0* cij_kl(8)* costhetasq + &
- (cij_kl(9) - cij_kl(20))* sinphi* sintwotheta))
-
- cij_kll(9) = cij_kl(11)* cosphifour* sintheta - sinphi*sinphisq* (cij_kl(5)* costheta + &
- cij_kl(6)* sinphi* sintheta) + cosphisq* sinphi* (-(cij_kl(10) + &
- cij_kl(18))* costheta + &
- 3.d0* (cij_kl(6) - cij_kl(11))* sinphi* sintheta) + &
- cosphi* sinphisq* ((cij_kl(4) + cij_kl(20))* costheta + &
- 2.d0* (-2.d0* cij_kl(1) + cij_kl(2) + cij_kl(21))* sinphi* sintheta) + &
- cosphi*cosphisq* (cij_kl(9)* costheta - 2.d0* (cij_kl(2) - 2.d0* cij_kl(7) + &
- cij_kl(21))* sinphi* sintheta)
-
- cij_kll(10) = 1.d0/4.d0* (4.d0* costwotheta* (cij_kl(10)* cosphi*cosphisq + &
- (cij_kl(9) - cij_kl(20))* cosphisq* sinphi + &
- (cij_kl(5) - cij_kl(18))* cosphi* sinphisq + &
- cij_kl(4)* sinphi*sinphisq) + (cij_kl(1) + 3.d0* cij_kl(2) - &
- 2.d0* cij_kl(3) + cij_kl(7) - &
- 2.d0* cij_kl(8) - cij_kl(21) + 2.d0* (cij_kl(3) - cij_kl(8))* costwophi + &
- (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* cosfourphi + &
- 2.d0* cij_kl(15)* sintwophi + &
- (-cij_kl(6) + cij_kl(11))* sinfourphi)* sintwotheta)
-
- cij_kll(11) = 1.d0/4.d0* (2.d0* costheta* ((cij_kl(6) + cij_kl(11))* costwophi + &
- (-cij_kl(6) + cij_kl(11))* cosfourphi + &
- 2.d0* (-cij_kl(1) + cij_kl(7))* sintwophi + &
- (cij_kl(1) - cij_kl(2) + cij_kl(7) - cij_kl(21))* sinfourphi) + &
- (-(cij_kl(4) + 3.d0* cij_kl(9) + cij_kl(20))* cosphi + &
- (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi + &
- (3.d0* cij_kl(5) + cij_kl(10) + cij_kl(18))* sinphi + &
- (-cij_kl(5) + cij_kl(10) + cij_kl(18))* sinthreephi)* sintheta)
-
- cij_kll(12) = 1.d0/16.d0* (cij_kl(16) - 2.d0* cij_kl(16)* cosfourtheta* sinphisq + &
- costwophi* (-cij_kl(16) + 8.d0* costheta* sinthetasq* ((cij_kl(3) - &
- cij_kl(8) + cij_kl(19))* costheta + &
- (cij_kl(5) - cij_kl(10) - cij_kl(18))* cosphi* sintheta)) + &
- 2.d0* (cij_kl(15) + cij_kl(17))* sintwophi* sintwothetasq + &
- 2.d0* (8.d0* cij_kl(12)* costhetafour + &
- 8.d0* cij_kl(14)* cosphi* costheta*costhetasq* sintheta + &
- 4.d0* cosphi* costheta* (cij_kl(5) + cij_kl(10) + cij_kl(18) + &
- (cij_kl(4) + cij_kl(20))* sintwophi)* &
- sintheta*sinthetasq + 8.d0* cij_kl(1)* cosphifour* sinthetafour + &
- 8.d0* cij_kl(6)* cosphi*cosphisq* sinphi* sinthetafour + &
- 8.d0* cij_kl(11)* cosphi* sinphi*sinphisq* sinthetafour + &
- 8.d0* cij_kl(7)* sinphifour* sinthetafour + &
- 2.d0* cij_kl(2)* sintwophisq* sinthetafour + &
- 2.d0* cij_kl(21)* sintwophisq* sinthetafour + &
- 2.d0* cij_kl(13)* sinphi* sintwotheta + &
- 2.d0* cij_kl(9)* sinphi*sinphisq* sintwotheta + &
- cij_kl(3)* sintwothetasq + cij_kl(8)* sintwothetasq + &
- cij_kl(19)* sintwothetasq + cij_kl(13)* sinphi* sinfourtheta - &
- cij_kl(9)* sinphi*sinphisq* sinfourtheta))
-
- cij_kll(13) = 1.d0/8.d0* (cosphi* costheta* (cij_kl(4) + 3.d0* cij_kl(9) + &
- 4.d0* cij_kl(13) + cij_kl(20) - (cij_kl(4) + 3.d0* cij_kl(9) - &
- 4.d0* cij_kl(13) + cij_kl(20))* costwotheta) + 4.d0* (-cij_kl(1) - &
- cij_kl(3) + cij_kl(7) + cij_kl(8) + cij_kl(16) - cij_kl(19) + &
- (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + cij_kl(16) - &
- cij_kl(19))* costwotheta)* sintwophi* sintheta + &
- 4.d0* (cij_kl(6) - cij_kl(11))* cosfourphi* sinthetasq*sintheta - &
- 4.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
- cij_kl(21))* sinfourphi* sinthetasq*sintheta + &
- costheta* ((-3.d0* cij_kl(5) - cij_kl(10) - 4.d0* cij_kl(14) - &
- cij_kl(18) + (3.d0* cij_kl(5) + cij_kl(10) - 4.d0* cij_kl(14) + &
- cij_kl(18))* costwotheta)* sinphi + 6.d0* ((cij_kl(4) - cij_kl(9) + &
- cij_kl(20))* costhreephi + (-cij_kl(5) + cij_kl(10) + &
- cij_kl(18))* sinthreephi)* sinthetasq) + costwophi* ((3* cij_kl(6) + &
- 3.d0* cij_kl(11) + 2.d0* (cij_kl(15) + cij_kl(17)))* sintheta - &
- (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + &
- cij_kl(17)))* sinthreetheta))
-
- cij_kll(14) = 1.d0/4.d0* (2.d0* cij_kl(13)* (costwotheta + cosfourtheta)* sinphi + &
- 8.d0* costheta*costhetasq* (-2.d0* cij_kl(12) + cij_kl(8)* sinphisq)* sintheta + &
- 4.d0* (cij_kl(4) + cij_kl(20))* cosphisq* (1.d0 + &
- 2.d0* costwotheta)* sinphi* sinthetasq + &
- 4.d0* cij_kl(9)* (1.d0 + 2.d0* costwotheta)* sinphi*sinphisq* sinthetasq + &
- 16.d0* cij_kl(1)* cosphifour* costheta* sintheta*sinthetasq + &
- 4.d0* costheta* (-2.d0* cij_kl(8)* sinphisq + 4.d0* cij_kl(7)* sinphifour + &
- (cij_kl(2) + cij_kl(21))* sintwophisq)* sintheta*sinthetasq + &
- 4.d0* cosphi*cosphisq* sinthetasq* (cij_kl(5) + 2.d0* cij_kl(5)* costwotheta + &
- 4.d0* cij_kl(6)* costheta* sinphi* sintheta) + &
- 2.d0* cosphi* (cosfourtheta* (cij_kl(14) - (cij_kl(10) + cij_kl(18))* sinphisq) + &
- costwotheta* (cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq) + &
- 8.d0* cij_kl(11)* costheta* sinphi*sinphisq* sintheta*sinthetasq) + &
- (cij_kl(3) + cij_kl(16) + cij_kl(19) + (cij_kl(3) - cij_kl(16) + &
- cij_kl(19))* costwophi + (cij_kl(15) + cij_kl(17))* sintwophi)* sinfourtheta)
-
- cij_kll(15) = costwophi* costheta* (-cij_kl(17) + (cij_kl(15) + cij_kl(17))* costhetasq) + &
- 1.d0/16.d0* (-((11.d0* cij_kl(4) + cij_kl(9) + 4.d0* cij_kl(13) - &
- 5.d0* cij_kl(20))* cosphi + (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi - &
- (cij_kl(5) + 11.d0* cij_kl(10) + 4.d0* cij_kl(14) - &
- 5.d0* cij_kl(18))* sinphi + (-cij_kl(5) + cij_kl(10) + &
- cij_kl(18))* sinthreephi)* sintheta + &
- 8.d0* costheta* ((-cij_kl(1) - cij_kl(3) + cij_kl(7) + cij_kl(8) - cij_kl(16) +&
- cij_kl(19) + (cij_kl(1) - cij_kl(3) - &
- cij_kl(7) + cij_kl(8) + cij_kl(16) - cij_kl(19))* costwotheta)* sintwophi +&
- ((cij_kl(6) + cij_kl(11))* costwophi + &
- (cij_kl(6) - cij_kl(11))* cosfourphi + (-cij_kl(1) + cij_kl(2) - cij_kl(7) +&
- cij_kl(21))* sinfourphi)* sinthetasq) +&
- ((cij_kl(4) + 3.d0* cij_kl(9) - 4.d0* cij_kl(13) + cij_kl(20))* cosphi + &
- 3.d0* (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi - &
- (3.d0* cij_kl(5) + cij_kl(10) - 4.d0* cij_kl(14) + cij_kl(18))* sinphi + &
- 3.d0* (-cij_kl(5) + cij_kl(10) + cij_kl(18))* sinthreephi)* sinthreetheta)
-
- cij_kll(16) = 1.d0/4.d0*(cij_kl(1) - cij_kl(2) + cij_kl(7) + cij_kl(16) + &
- cij_kl(19) + cij_kl(21) + 2.d0*(cij_kl(16) - cij_kl(19))*costwophi* costhetasq + &
- (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(16) + &
- cij_kl(19) - cij_kl(21))*costwotheta - 2.d0* cij_kl(17)* costhetasq* sintwophi + &
- 2.d0* ((-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* cosfourphi + &
- (-cij_kl(6) + cij_kl(11))* sinfourphi)* sinthetasq + ((cij_kl(5) - cij_kl(10) +&
- cij_kl(18))* cosphi + (-cij_kl(5) + cij_kl(10) + cij_kl(18))* costhreephi +&
- (-cij_kl(4) + cij_kl(9) + cij_kl(20))* sinphi - &
- (cij_kl(4) - cij_kl(9) + cij_kl(20))* sinthreephi)* sintwotheta)
-
- cij_kll(17) = 1.d0/8.d0* (4.d0* costwophi* costheta* (cij_kl(6) + cij_kl(11) - &
- 2.d0* cij_kl(15) - (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + &
- cij_kl(17)))* costwotheta) - (2.d0* cosphi* (-3.d0* cij_kl(4) +&
- cij_kl(9) + 2.d0* cij_kl(13) + cij_kl(20) + (cij_kl(4) - cij_kl(9) + &
- cij_kl(20))* costwophi) - (cij_kl(5) - 5.d0* cij_kl(10) + &
- 4.d0* cij_kl(14) + 3.d0* cij_kl(18))* sinphi + (-cij_kl(5) + cij_kl(10) + &
- cij_kl(18))* sinthreephi)* sintheta + &
- 8.d0* costheta* ((-cij_kl(1) + cij_kl(3) + cij_kl(7) - cij_kl(8) + &
- (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + cij_kl(16) - &
- cij_kl(19))* costwotheta)* sintwophi + ((cij_kl(6) - cij_kl(11))* cosfourphi + &
- (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* sinfourphi)* sinthetasq) +&
- ((cij_kl(4) + 3.d0* cij_kl(9) - 4.d0* cij_kl(13) + cij_kl(20))* cosphi + &
- 3.d0* (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi - &
- (3.d0* cij_kl(5) + cij_kl(10) - 4.d0* cij_kl(14) + cij_kl(18))* sinphi + &
- 3.d0* (-cij_kl(5) + cij_kl(10) + cij_kl(18))* sinthreephi)* sinthreetheta)
-
- cij_kll(18) = 1.d0/2.d0* ((cij_kl(5) - cij_kl(10) + cij_kl(18))* cosphi* costwotheta - &
- (cij_kl(5) - cij_kl(10) - cij_kl(18))* costhreephi* costwotheta - &
- 2.d0* (cij_kl(4) - cij_kl(9) + &
- (cij_kl(4) - cij_kl(9) + cij_kl(20))* costwophi)* costwotheta* sinphi + &
- (cij_kl(1) - cij_kl(2) + cij_kl(7) - cij_kl(16) - cij_kl(19) + cij_kl(21) + &
- (-cij_kl(16) + cij_kl(19))* costwophi + &
- (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* cosfourphi + &
- cij_kl(17)* sintwophi + &
- (-cij_kl(6) + cij_kl(11))* sinfourphi)* sintwotheta)
-
- cij_kll(19) = 1.d0/4.d0* (cij_kl(16) - cij_kl(16)* costwophi + &
- (-cij_kl(15) + cij_kl(17))* sintwophi + &
- 4.d0* cij_kl(12)* sintwothetasq + &
- 2.d0* (2.d0* cij_kl(1)* cosphifour* sintwothetasq + &
- cosphi*cosphisq* (8.d0* cij_kl(6)* costhetasq* sinphi* sinthetasq + &
- cij_kl(5)* sinfourtheta) + cosphisq* (-cij_kl(3) + cij_kl(19) + (cij_kl(3) +&
- cij_kl(19))* cosfourtheta + (cij_kl(4) + cij_kl(20))* sinphi* sinfourtheta) + &
- sinphi* (cosfourtheta* ((cij_kl(15) + cij_kl(17))* cosphi + &
- cij_kl(16)* sinphi) + (cij_kl(2) + cij_kl(7) - 2.d0* cij_kl(8) + cij_kl(21) + &
- (cij_kl(2) - cij_kl(7) + cij_kl(21))* costwophi)* sinphi* sintwothetasq + &
- (-cij_kl(13) + cij_kl(9)* sinphisq)* sinfourtheta) + &
- cosphi* (8.d0* cij_kl(11)* costhetasq* sinphi*sinphisq* sinthetasq + &
- (-cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq)* sinfourtheta)))
-
- cij_kll(20) = 1.d0/8.d0* (2.d0* cosphi* costheta* (-3.d0* cij_kl(4) - cij_kl(9) + &
- 4.d0* cij_kl(13) + cij_kl(20) + (cij_kl(4) + 3.d0* cij_kl(9) - &
- 4.d0* cij_kl(13) + cij_kl(20))* costwotheta) + &
- (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi* (costheta + &
- 3.d0* costhreetheta) - &
- 2.d0* costheta* (-cij_kl(5) - 3.d0* cij_kl(10) + 4.d0* cij_kl(14) + &
- cij_kl(18) + (3.d0* cij_kl(5) + &
- cij_kl(10) - 4.d0* cij_kl(14) + cij_kl(18))*costwotheta)* sinphi - &
- (cij_kl(5) - cij_kl(10) - cij_kl(18))* &
- (costheta + 3.d0* costhreetheta)* sinthreephi + 8.d0* (cij_kl(6) - &
- cij_kl(11))* cosfourphi* costhetasq* sintheta - 8.d0* (cij_kl(1) - &
- cij_kl(3) - cij_kl(7) + cij_kl(8) + &
- (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + cij_kl(16) - &
- cij_kl(19))* costwotheta)* sintwophi* sintheta - &
- 8.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
- cij_kl(21))* costhetasq* sinfourphi* sintheta + &
- 2.d0* costwophi* ((cij_kl(6) + cij_kl(11) - 2.d0* cij_kl(15) + &
- 2.d0* cij_kl(17))* sintheta + &
- (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + cij_kl(17)))* sinthreetheta))
-
- cij_kll(21) = 1.d0/4.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) + cij_kl(16) + &
- cij_kl(19) + cij_kl(21) - 2.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
- cij_kl(21))* cosfourphi* costhetasq + &
- (cij_kl(1) - cij_kl(2) + cij_kl(7) - cij_kl(16) - cij_kl(19) + &
- cij_kl(21))* costwotheta + &
- 2.d0* (-cij_kl(6) + cij_kl(11))* costhetasq* sinfourphi - &
- 2.d0* ((-cij_kl(16) + cij_kl(19))* costwophi + cij_kl(17)* sintwophi)* sinthetasq - &
- ((cij_kl(5) - cij_kl(10) + cij_kl(18))* cosphi + (-cij_kl(5) + cij_kl(10) +&
- cij_kl(18))* costhreephi + &
- (-cij_kl(4) + cij_kl(9) + cij_kl(20))* sinphi - (cij_kl(4) - cij_kl(9) + &
- cij_kl(20))* sinthreephi)* sintwotheta)
-
- end subroutine rotate_kernels_dble
-
-!-----------------------------------------------------------------------------
-
subroutine compute_kernels_hessian(ibool_crust_mantle, &
hess_kl_crust_mantle, &
accel_crust_mantle,b_accel_crust_mantle, &
@@ -870,25 +519,41 @@
real(kind=CUSTOM_REAL) deltat
! local parameters
- integer :: i,j,k,ispec,iglob
+ integer :: ispec,iglob
- ! crust_mantle
- do ispec = 1, NSPEC_CRUST_MANTLE
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool_crust_mantle(i,j,k,ispec)
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#else
+ integer :: i,j,k
+#endif
- ! approximates hessian
- ! term with adjoint acceleration and backward/reconstructed acceleration
- hess_kl_crust_mantle(i,j,k,ispec) = hess_kl_crust_mantle(i,j,k,ispec) &
- + deltat * (accel_crust_mantle(1,iglob) * b_accel_crust_mantle(1,iglob) &
- + accel_crust_mantle(2,iglob) * b_accel_crust_mantle(2,iglob) &
- + accel_crust_mantle(3,iglob) * b_accel_crust_mantle(3,iglob) )
+! approximates Hessian term with adjoint acceleration and backward/reconstructed acceleration
+ ! crust_mantle
+ do ispec = 1, NSPEC_CRUST_MANTLE
+
+#ifdef FORCE_VECTORIZATION
+ do ijk = 1,NGLLCUBE
+ iglob = ibool_crust_mantle(ijk,1,1,ispec)
+ hess_kl_crust_mantle(ijk,1,1,ispec) = hess_kl_crust_mantle(ijk,1,1,ispec) &
+ + deltat * (accel_crust_mantle(1,iglob) * b_accel_crust_mantle(1,iglob) &
+ + accel_crust_mantle(2,iglob) * b_accel_crust_mantle(2,iglob) &
+ + accel_crust_mantle(3,iglob) * b_accel_crust_mantle(3,iglob))
enddo
+#else
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ hess_kl_crust_mantle(i,j,k,ispec) = hess_kl_crust_mantle(i,j,k,ispec) &
+ + deltat * (accel_crust_mantle(1,iglob) * b_accel_crust_mantle(1,iglob) &
+ + accel_crust_mantle(2,iglob) * b_accel_crust_mantle(2,iglob) &
+ + accel_crust_mantle(3,iglob) * b_accel_crust_mantle(3,iglob))
+ enddo
enddo
enddo
+#endif
+
enddo
end subroutine compute_kernels_hessian
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk 2013-07-13 14:03:38 UTC (rev 22583)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk 2013-07-13 14:08:06 UTC (rev 22584)
@@ -183,8 +183,8 @@
$O/compute_forces_inner_core_Dev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_inner_core_Dev.F90
${FCCOMPILE_CHECK} -c -o $O/compute_forces_inner_core_Dev.o ${FCFLAGS_f90} $S/compute_forces_inner_core_Dev.F90
-$O/compute_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_kernels.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o ${FCFLAGS_f90} $S/compute_kernels.f90
+$O/compute_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_kernels.F90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o ${FCFLAGS_f90} $S/compute_kernels.F90
$O/compute_seismograms.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_seismograms.f90
${FCCOMPILE_CHECK} -c -o $O/compute_seismograms.o ${FCFLAGS_f90} $S/compute_seismograms.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_kernels.f90 2013-07-13 14:03:38 UTC (rev 22583)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_kernels.f90 2013-07-13 14:08:06 UTC (rev 22584)
@@ -859,3 +859,381 @@
end subroutine save_kernels_hessian
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine rotate_kernels_dble(cij_kl,cij_kll,theta_in,phi_in)
+
+! Purpose : compute the kernels in r,theta,phi (cij_kll)
+! from the kernels in x,y,z (cij_kl) (x,y,z <-> r,theta,phi)
+! At r,theta,phi fixed
+! theta and phi are in radians
+
+! Coeff from Min's routine rotate_anisotropic_tensor
+! with the help of Collect[Expand[cij],{dij}] in Mathematica
+
+! Definition of the output array cij_kll :
+! cij_kll(1) = C11 ; cij_kll(2) = C12 ; cij_kll(3) = C13
+! cij_kll(4) = C14 ; cij_kll(5) = C15 ; cij_kll(6) = C16
+! cij_kll(7) = C22 ; cij_kll(8) = C23 ; cij_kll(9) = C24
+! cij_kll(10) = C25 ; cij_kll(11) = C26 ; cij_kll(12) = C33
+! cij_kll(13) = C34 ; cij_kll(14) = C35 ; cij_kll(15) = C36
+! cij_kll(16) = C44 ; cij_kll(17) = C45 ; cij_kll(18) = C46
+! cij_kll(19) = C55 ; cij_kll(20) = C56 ; cij_kll(21) = C66
+! where the Cij (Voigt's notation) are defined as function of
+! the components of the elastic tensor in spherical coordinates
+! by eq. (A.1) of Chen & Tromp, GJI 168 (2007)
+
+ implicit none
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL), intent(in) :: theta_in,phi_in
+
+ real(kind=CUSTOM_REAL), dimension(21), intent(in) :: cij_kl
+ real(kind=CUSTOM_REAL), dimension(21), intent(out) :: cij_kll
+
+ double precision :: theta,phi
+ double precision :: costheta,sintheta,cosphi,sinphi
+ double precision :: costhetasq,sinthetasq,cosphisq,sinphisq
+ double precision :: costwotheta,sintwotheta,costwophi,sintwophi
+ double precision :: cosfourtheta,sinfourtheta,cosfourphi,sinfourphi
+ double precision :: costhetafour,sinthetafour,cosphifour,sinphifour
+ double precision :: sintwophisq,sintwothetasq
+ double precision :: costhreetheta,sinthreetheta,costhreephi,sinthreephi
+
+ if (CUSTOM_REAL == SIZE_REAL) then
+ theta = dble(theta_in)
+ phi = dble(phi_in)
+ else
+ theta = theta_in
+ phi = phi_in
+ endif
+
+ costheta = dcos(theta)
+ sintheta = dsin(theta)
+ cosphi = dcos(phi)
+ sinphi = dsin(phi)
+
+ costhetasq = costheta * costheta
+ sinthetasq = sintheta * sintheta
+ cosphisq = cosphi * cosphi
+ sinphisq = sinphi * sinphi
+
+ costhetafour = costhetasq * costhetasq
+ sinthetafour = sinthetasq * sinthetasq
+ cosphifour = cosphisq * cosphisq
+ sinphifour = sinphisq * sinphisq
+
+ costwotheta = dcos(2.d0*theta)
+ sintwotheta = dsin(2.d0*theta)
+ costwophi = dcos(2.d0*phi)
+ sintwophi = dsin(2.d0*phi)
+
+ costhreetheta=dcos(3.d0*theta)
+ sinthreetheta=dsin(3.d0*theta)
+ costhreephi=dcos(3.d0*phi)
+ sinthreephi=dsin(3.d0*phi)
+
+ cosfourtheta = dcos(4.d0*theta)
+ sinfourtheta = dsin(4.d0*theta)
+ cosfourphi = dcos(4.d0*phi)
+ sinfourphi = dsin(4.d0*phi)
+ sintwothetasq = sintwotheta * sintwotheta
+ sintwophisq = sintwophi * sintwophi
+
+ cij_kll(1) = ONE_SIXTEENTH* (cij_kl(16) - cij_kl(16)* costwophi + &
+ 16.d0* cosphi*cosphisq* costhetafour* (cij_kl(1)* cosphi + cij_kl(6)* sinphi) + &
+ 2.d0* (cij_kl(15) + cij_kl(17))* sintwophi* sintwothetasq - &
+ 2.d0* (cij_kl(16)* cosfourtheta* sinphisq + &
+ 2.d0* costhetafour* (-4* cij_kl(7)* sinphifour - &
+ (cij_kl(2) + cij_kl(21))* sintwophisq) + &
+ 8.d0* cij_kl(5)* cosphi*cosphisq* costheta*costhetasq* sintheta - &
+ 8.d0* cij_kl(8)* costhetasq* sinphisq* sinthetasq - &
+ 8.d0* cij_kl(12)* sinthetafour + &
+ 8.d0* cosphisq* costhetasq* sintheta* ((cij_kl(4) + &
+ cij_kl(20))* costheta* sinphi - &
+ (cij_kl(3) + cij_kl(19))*sintheta) + &
+ 8.d0* cosphi* costheta* (-cij_kl(11)* costheta*costhetasq* &
+ sinphi*sinphisq + (cij_kl(10) + cij_kl(18))* costhetasq* sinphisq* sintheta + &
+ cij_kl(14)* sintheta*sinthetasq) + 2.d0* sinphi* (cij_kl(13) + &
+ cij_kl(9)* sinphisq)* sintwotheta + &
+ sinphi* (-cij_kl(13) + cij_kl(9)* sinphisq)* sinfourtheta))
+
+ cij_kll(2) = ONE_FOURTH* (costhetasq* (cij_kl(1) + 3.d0* cij_kl(2) + cij_kl(7) - &
+ cij_kl(21) + (-cij_kl(1) + cij_kl(2) - cij_kl(7) + &
+ cij_kl(21))* cosfourphi + (-cij_kl(6) + cij_kl(11))* sinfourphi) + &
+ 4.d0* (cij_kl(8)* cosphisq - cij_kl(15)* cosphi* sinphi + &
+ cij_kl(3)* sinphisq)* sinthetasq - &
+ 2.d0* (cij_kl(10)* cosphisq*cosphi + &
+ (cij_kl(9) - cij_kl(20))* cosphisq* sinphi + &
+ (cij_kl(5) - cij_kl(18))* cosphi* sinphisq + &
+ cij_kl(4)* sinphisq*sinphi)* sintwotheta)
+
+ cij_kll(3) = ONE_EIGHTH* (sintwophi* (3.d0* cij_kl(15) - cij_kl(17) + &
+ 4.d0* (cij_kl(2) + cij_kl(21))* costhetasq* sintwophi* sinthetasq) + &
+ 4.d0* cij_kl(12)* sintwothetasq + 4.d0* cij_kl(1)* cosphifour* sintwothetasq + &
+ 2.d0* cosphi*cosphisq* (8.d0* cij_kl(6)* costhetasq* sinphi* sinthetasq + &
+ cij_kl(5)* sinfourtheta) + 2.d0* cosphisq* (3.d0* cij_kl(3) - cij_kl(19) + &
+ (cij_kl(3) + cij_kl(19))* cosfourtheta + &
+ (cij_kl(4) + cij_kl(20))* sinphi* sinfourtheta) + &
+ 2.d0* sinphi* (sinphi* (3.d0* cij_kl(8) - &
+ cij_kl(16) + (cij_kl(8) + cij_kl(16))* cosfourtheta + &
+ 2.d0* cij_kl(7)* sinphisq* sintwothetasq)+ &
+ (-cij_kl(13) + cij_kl(9)* sinphisq)* sinfourtheta)+ &
+ 2.d0* cosphi* ((cij_kl(15) + cij_kl(17))* cosfourtheta* sinphi + &
+ 8.d0* cij_kl(11)* costhetasq* sinphi*sinphisq* sinthetasq + &
+ (-cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq)*sinfourtheta))
+
+ cij_kll(4) = ONE_EIGHTH* (cosphi* costheta *(5.d0* cij_kl(4) - &
+ cij_kl(9) + 4.d0* cij_kl(13) - &
+ 3.d0* cij_kl(20) + (cij_kl(4) + 3.d0* cij_kl(9) - &
+ 4.d0* cij_kl(13) + cij_kl(20))* costwotheta) + &
+ ONE_HALF* (cij_kl(4) - cij_kl(9) + &
+ cij_kl(20))* costhreephi * (costheta + 3.d0* costhreetheta) - &
+ costheta* (-cij_kl(5) + 5.d0* cij_kl(10) + &
+ 4.d0* cij_kl(14) - 3.d0* cij_kl(18) + &
+ (3.d0* cij_kl(5) + cij_kl(10) - &
+ 4.d0* cij_kl(14) + cij_kl(18))* costwotheta)* sinphi - &
+ ONE_HALF* (cij_kl(5) - cij_kl(10) - cij_kl(18))* (costheta + &
+ 3.d0* costhreetheta)* sinthreephi + &
+ 4.d0* (cij_kl(6) - cij_kl(11))* cosfourphi* costhetasq* sintheta - &
+ 4.d0* (cij_kl(1) + cij_kl(3) - cij_kl(7) - cij_kl(8) + cij_kl(16) - cij_kl(19) + &
+ (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + &
+ cij_kl(16) - cij_kl(19))* costwotheta)* sintwophi* sintheta - &
+ 4.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
+ cij_kl(21))* costhetasq* sinfourphi* sintheta + &
+ costwophi* ((cij_kl(6) + cij_kl(11) + 6.d0* cij_kl(15) - &
+ 2.d0* cij_kl(17))* sintheta + &
+ (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + cij_kl(17)))* sinthreetheta))
+
+ cij_kll(5) = ONE_FOURTH* (2.d0* (cij_kl(4) + &
+ cij_kl(20))* cosphisq* (costwotheta + cosfourtheta)* sinphi + &
+ 2.d0* cij_kl(9)* (costwotheta + cosfourtheta)* sinphi*sinphisq + &
+ 16.d0* cij_kl(1)* cosphifour* costheta*costhetasq* sintheta + &
+ 4.d0* costheta*costhetasq* (-2.d0* cij_kl(8)* sinphisq + &
+ 4.d0* cij_kl(7)* sinphifour + &
+ (cij_kl(2) + cij_kl(21))* sintwophisq)* sintheta + &
+ 4.d0* cij_kl(13)* (1.d0 + 2.d0* costwotheta)* sinphi* sinthetasq + &
+ 8.d0* costheta* (-2.d0* cij_kl(12) + cij_kl(8)* sinphisq)* sintheta*sinthetasq + &
+ 2.d0* cosphi*cosphisq* (cij_kl(5)* (costwotheta + cosfourtheta) + &
+ 8.d0* cij_kl(6)* costheta*costhetasq* sinphi* sintheta) + &
+ 2.d0* cosphi* (cosfourtheta* (-cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq) + &
+ costwotheta* (cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq) + &
+ 8.d0* cij_kl(11)* costheta*costhetasq* sinphi*sinphisq* sintheta) - &
+ (cij_kl(3) + cij_kl(16) + cij_kl(19) + &
+ (cij_kl(3) - cij_kl(16) + cij_kl(19))* costwophi + &
+ (cij_kl(15) + cij_kl(17))* sintwophi)* sinfourtheta)
+
+ cij_kll(6) = ONE_HALF* costheta*costhetasq* ((cij_kl(6) + cij_kl(11))* costwophi + &
+ (cij_kl(6) - cij_kl(11))* cosfourphi + 2.d0* (-cij_kl(1) + cij_kl(7))* sintwophi + &
+ (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* sinfourphi) + &
+ ONE_FOURTH* costhetasq* (-(cij_kl(4) + 3* cij_kl(9) + cij_kl(20))* cosphi - &
+ 3.d0* (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi + &
+ (3.d0* cij_kl(5) + cij_kl(10) + cij_kl(18))* sinphi + &
+ 3.d0* (cij_kl(5) - cij_kl(10) - cij_kl(18))* sinthreephi)* sintheta + &
+ costheta* ((cij_kl(15) + cij_kl(17))* costwophi + &
+ (-cij_kl(3) + cij_kl(8) + cij_kl(16) - cij_kl(19))* sintwophi)* sinthetasq + &
+ (-cij_kl(13)* cosphi + cij_kl(14)* sinphi)* sintheta*sinthetasq
+
+ cij_kll(7) = cij_kl(7)* cosphifour - cij_kl(11)* cosphi*cosphisq* sinphi + &
+ (cij_kl(2) + cij_kl(21))* cosphisq* sinphisq - &
+ cij_kl(6)* cosphi* sinphi*sinphisq + &
+ cij_kl(1)* sinphifour
+
+ cij_kll(8) = ONE_HALF* (2.d0* costhetasq* sinphi* (-cij_kl(15)* cosphi + &
+ cij_kl(3)* sinphi) + 2.d0* cij_kl(2)* cosphifour* sinthetasq + &
+ (2.d0* cij_kl(2)* sinphifour + &
+ (cij_kl(1) + cij_kl(7) - cij_kl(21))* sintwophisq)* sinthetasq + &
+ cij_kl(4)* sinphi*sinphisq* sintwotheta + &
+ cosphi*cosphisq* (2.d0* (-cij_kl(6) + cij_kl(11))* sinphi* sinthetasq + &
+ cij_kl(10)* sintwotheta) + cosphi* sinphisq* (2.d0* (cij_kl(6) - &
+ cij_kl(11))* sinphi* sinthetasq + &
+ (cij_kl(5) - cij_kl(18))* sintwotheta) + &
+ cosphisq* (2.d0* cij_kl(8)* costhetasq + &
+ (cij_kl(9) - cij_kl(20))* sinphi* sintwotheta))
+
+ cij_kll(9) = cij_kl(11)* cosphifour* sintheta - sinphi*sinphisq* (cij_kl(5)* costheta + &
+ cij_kl(6)* sinphi* sintheta) + cosphisq* sinphi* (-(cij_kl(10) + &
+ cij_kl(18))* costheta + &
+ 3.d0* (cij_kl(6) - cij_kl(11))* sinphi* sintheta) + &
+ cosphi* sinphisq* ((cij_kl(4) + cij_kl(20))* costheta + &
+ 2.d0* (-2.d0* cij_kl(1) + cij_kl(2) + cij_kl(21))* sinphi* sintheta) + &
+ cosphi*cosphisq* (cij_kl(9)* costheta - 2.d0* (cij_kl(2) - 2.d0* cij_kl(7) + &
+ cij_kl(21))* sinphi* sintheta)
+
+ cij_kll(10) = ONE_FOURTH* (4.d0* costwotheta* (cij_kl(10)* cosphi*cosphisq + &
+ (cij_kl(9) - cij_kl(20))* cosphisq* sinphi + &
+ (cij_kl(5) - cij_kl(18))* cosphi* sinphisq + &
+ cij_kl(4)* sinphi*sinphisq) + (cij_kl(1) + 3.d0* cij_kl(2) - &
+ 2.d0* cij_kl(3) + cij_kl(7) - &
+ 2.d0* cij_kl(8) - cij_kl(21) + 2.d0* (cij_kl(3) - cij_kl(8))* costwophi + &
+ (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* cosfourphi + &
+ 2.d0* cij_kl(15)* sintwophi + &
+ (-cij_kl(6) + cij_kl(11))* sinfourphi)* sintwotheta)
+
+ cij_kll(11) = ONE_FOURTH* (2.d0* costheta* ((cij_kl(6) + cij_kl(11))* costwophi + &
+ (-cij_kl(6) + cij_kl(11))* cosfourphi + &
+ 2.d0* (-cij_kl(1) + cij_kl(7))* sintwophi + &
+ (cij_kl(1) - cij_kl(2) + cij_kl(7) - cij_kl(21))* sinfourphi) + &
+ (-(cij_kl(4) + 3.d0* cij_kl(9) + cij_kl(20))* cosphi + &
+ (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi + &
+ (3.d0* cij_kl(5) + cij_kl(10) + cij_kl(18))* sinphi + &
+ (-cij_kl(5) + cij_kl(10) + cij_kl(18))* sinthreephi)* sintheta)
+
+ cij_kll(12) = ONE_SIXTEENTH* (cij_kl(16) - 2.d0* cij_kl(16)* cosfourtheta* sinphisq + &
+ costwophi* (-cij_kl(16) + 8.d0* costheta* sinthetasq* ((cij_kl(3) - &
+ cij_kl(8) + cij_kl(19))* costheta + &
+ (cij_kl(5) - cij_kl(10) - cij_kl(18))* cosphi* sintheta)) + &
+ 2.d0* (cij_kl(15) + cij_kl(17))* sintwophi* sintwothetasq + &
+ 2.d0* (8.d0* cij_kl(12)* costhetafour + &
+ 8.d0* cij_kl(14)* cosphi* costheta*costhetasq* sintheta + &
+ 4.d0* cosphi* costheta* (cij_kl(5) + cij_kl(10) + cij_kl(18) + &
+ (cij_kl(4) + cij_kl(20))* sintwophi)* &
+ sintheta*sinthetasq + 8.d0* cij_kl(1)* cosphifour* sinthetafour + &
+ 8.d0* cij_kl(6)* cosphi*cosphisq* sinphi* sinthetafour + &
+ 8.d0* cij_kl(11)* cosphi* sinphi*sinphisq* sinthetafour + &
+ 8.d0* cij_kl(7)* sinphifour* sinthetafour + &
+ 2.d0* cij_kl(2)* sintwophisq* sinthetafour + &
+ 2.d0* cij_kl(21)* sintwophisq* sinthetafour + &
+ 2.d0* cij_kl(13)* sinphi* sintwotheta + &
+ 2.d0* cij_kl(9)* sinphi*sinphisq* sintwotheta + &
+ cij_kl(3)* sintwothetasq + cij_kl(8)* sintwothetasq + &
+ cij_kl(19)* sintwothetasq + cij_kl(13)* sinphi* sinfourtheta - &
+ cij_kl(9)* sinphi*sinphisq* sinfourtheta))
+
+ cij_kll(13) = ONE_EIGHTH* (cosphi* costheta* (cij_kl(4) + 3.d0* cij_kl(9) + &
+ 4.d0* cij_kl(13) + cij_kl(20) - (cij_kl(4) + 3.d0* cij_kl(9) - &
+ 4.d0* cij_kl(13) + cij_kl(20))* costwotheta) + 4.d0* (-cij_kl(1) - &
+ cij_kl(3) + cij_kl(7) + cij_kl(8) + cij_kl(16) - cij_kl(19) + &
+ (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + cij_kl(16) - &
+ cij_kl(19))* costwotheta)* sintwophi* sintheta + &
+ 4.d0* (cij_kl(6) - cij_kl(11))* cosfourphi* sinthetasq*sintheta - &
+ 4.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
+ cij_kl(21))* sinfourphi* sinthetasq*sintheta + &
+ costheta* ((-3.d0* cij_kl(5) - cij_kl(10) - 4.d0* cij_kl(14) - &
+ cij_kl(18) + (3.d0* cij_kl(5) + cij_kl(10) - 4.d0* cij_kl(14) + &
+ cij_kl(18))* costwotheta)* sinphi + 6.d0* ((cij_kl(4) - cij_kl(9) + &
+ cij_kl(20))* costhreephi + (-cij_kl(5) + cij_kl(10) + &
+ cij_kl(18))* sinthreephi)* sinthetasq) + costwophi* ((3* cij_kl(6) + &
+ 3.d0* cij_kl(11) + 2.d0* (cij_kl(15) + cij_kl(17)))* sintheta - &
+ (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + &
+ cij_kl(17)))* sinthreetheta))
+
+ cij_kll(14) = ONE_FOURTH* (2.d0* cij_kl(13)* (costwotheta + cosfourtheta)* sinphi + &
+ 8.d0* costheta*costhetasq* (-2.d0* cij_kl(12) + cij_kl(8)* sinphisq)* sintheta + &
+ 4.d0* (cij_kl(4) + cij_kl(20))* cosphisq* (1.d0 + &
+ 2.d0* costwotheta)* sinphi* sinthetasq + &
+ 4.d0* cij_kl(9)* (1.d0 + 2.d0* costwotheta)* sinphi*sinphisq* sinthetasq + &
+ 16.d0* cij_kl(1)* cosphifour* costheta* sintheta*sinthetasq + &
+ 4.d0* costheta* (-2.d0* cij_kl(8)* sinphisq + 4.d0* cij_kl(7)* sinphifour + &
+ (cij_kl(2) + cij_kl(21))* sintwophisq)* sintheta*sinthetasq + &
+ 4.d0* cosphi*cosphisq* sinthetasq* (cij_kl(5) + 2.d0* cij_kl(5)* costwotheta + &
+ 4.d0* cij_kl(6)* costheta* sinphi* sintheta) + &
+ 2.d0* cosphi* (cosfourtheta* (cij_kl(14) - (cij_kl(10) + cij_kl(18))* sinphisq) + &
+ costwotheta* (cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq) + &
+ 8.d0* cij_kl(11)* costheta* sinphi*sinphisq* sintheta*sinthetasq) + &
+ (cij_kl(3) + cij_kl(16) + cij_kl(19) + (cij_kl(3) - cij_kl(16) + &
+ cij_kl(19))* costwophi + (cij_kl(15) + cij_kl(17))* sintwophi)* sinfourtheta)
+
+ cij_kll(15) = costwophi* costheta* (-cij_kl(17) + (cij_kl(15) + cij_kl(17))* costhetasq) + &
+ ONE_SIXTEENTH* (-((11.d0* cij_kl(4) + cij_kl(9) + 4.d0* cij_kl(13) - &
+ 5.d0* cij_kl(20))* cosphi + (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi - &
+ (cij_kl(5) + 11.d0* cij_kl(10) + 4.d0* cij_kl(14) - &
+ 5.d0* cij_kl(18))* sinphi + (-cij_kl(5) + cij_kl(10) + &
+ cij_kl(18))* sinthreephi)* sintheta + &
+ 8.d0* costheta* ((-cij_kl(1) - cij_kl(3) + cij_kl(7) + cij_kl(8) - cij_kl(16) +&
+ cij_kl(19) + (cij_kl(1) - cij_kl(3) - &
+ cij_kl(7) + cij_kl(8) + cij_kl(16) - cij_kl(19))* costwotheta)* sintwophi +&
+ ((cij_kl(6) + cij_kl(11))* costwophi + &
+ (cij_kl(6) - cij_kl(11))* cosfourphi + (-cij_kl(1) + cij_kl(2) - cij_kl(7) +&
+ cij_kl(21))* sinfourphi)* sinthetasq) +&
+ ((cij_kl(4) + 3.d0* cij_kl(9) - 4.d0* cij_kl(13) + cij_kl(20))* cosphi + &
+ 3.d0* (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi - &
+ (3.d0* cij_kl(5) + cij_kl(10) - 4.d0* cij_kl(14) + cij_kl(18))* sinphi + &
+ 3.d0* (-cij_kl(5) + cij_kl(10) + cij_kl(18))* sinthreephi)* sinthreetheta)
+
+ cij_kll(16) = ONE_FOURTH*(cij_kl(1) - cij_kl(2) + cij_kl(7) + cij_kl(16) + &
+ cij_kl(19) + cij_kl(21) + 2.d0*(cij_kl(16) - cij_kl(19))*costwophi* costhetasq + &
+ (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(16) + &
+ cij_kl(19) - cij_kl(21))*costwotheta - 2.d0* cij_kl(17)* costhetasq* sintwophi + &
+ 2.d0* ((-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* cosfourphi + &
+ (-cij_kl(6) + cij_kl(11))* sinfourphi)* sinthetasq + ((cij_kl(5) - cij_kl(10) +&
+ cij_kl(18))* cosphi + (-cij_kl(5) + cij_kl(10) + cij_kl(18))* costhreephi +&
+ (-cij_kl(4) + cij_kl(9) + cij_kl(20))* sinphi - &
+ (cij_kl(4) - cij_kl(9) + cij_kl(20))* sinthreephi)* sintwotheta)
+
+ cij_kll(17) = ONE_EIGHTH* (4.d0* costwophi* costheta* (cij_kl(6) + cij_kl(11) - &
+ 2.d0* cij_kl(15) - (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + &
+ cij_kl(17)))* costwotheta) - (2.d0* cosphi* (-3.d0* cij_kl(4) +&
+ cij_kl(9) + 2.d0* cij_kl(13) + cij_kl(20) + (cij_kl(4) - cij_kl(9) + &
+ cij_kl(20))* costwophi) - (cij_kl(5) - 5.d0* cij_kl(10) + &
+ 4.d0* cij_kl(14) + 3.d0* cij_kl(18))* sinphi + (-cij_kl(5) + cij_kl(10) + &
+ cij_kl(18))* sinthreephi)* sintheta + &
+ 8.d0* costheta* ((-cij_kl(1) + cij_kl(3) + cij_kl(7) - cij_kl(8) + &
+ (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + cij_kl(16) - &
+ cij_kl(19))* costwotheta)* sintwophi + ((cij_kl(6) - cij_kl(11))* cosfourphi + &
+ (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* sinfourphi)* sinthetasq) +&
+ ((cij_kl(4) + 3.d0* cij_kl(9) - 4.d0* cij_kl(13) + cij_kl(20))* cosphi + &
+ 3.d0* (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi - &
+ (3.d0* cij_kl(5) + cij_kl(10) - 4.d0* cij_kl(14) + cij_kl(18))* sinphi + &
+ 3.d0* (-cij_kl(5) + cij_kl(10) + cij_kl(18))* sinthreephi)* sinthreetheta)
+
+ cij_kll(18) = ONE_HALF* ((cij_kl(5) - cij_kl(10) + cij_kl(18))* cosphi* costwotheta - &
+ (cij_kl(5) - cij_kl(10) - cij_kl(18))* costhreephi* costwotheta - &
+ 2.d0* (cij_kl(4) - cij_kl(9) + &
+ (cij_kl(4) - cij_kl(9) + cij_kl(20))* costwophi)* costwotheta* sinphi + &
+ (cij_kl(1) - cij_kl(2) + cij_kl(7) - cij_kl(16) - cij_kl(19) + cij_kl(21) + &
+ (-cij_kl(16) + cij_kl(19))* costwophi + &
+ (-cij_kl(1) + cij_kl(2) - cij_kl(7) + cij_kl(21))* cosfourphi + &
+ cij_kl(17)* sintwophi + &
+ (-cij_kl(6) + cij_kl(11))* sinfourphi)* sintwotheta)
+
+ cij_kll(19) = ONE_FOURTH* (cij_kl(16) - cij_kl(16)* costwophi + &
+ (-cij_kl(15) + cij_kl(17))* sintwophi + &
+ 4.d0* cij_kl(12)* sintwothetasq + &
+ 2.d0* (2.d0* cij_kl(1)* cosphifour* sintwothetasq + &
+ cosphi*cosphisq* (8.d0* cij_kl(6)* costhetasq* sinphi* sinthetasq + &
+ cij_kl(5)* sinfourtheta) + cosphisq* (-cij_kl(3) + cij_kl(19) + (cij_kl(3) +&
+ cij_kl(19))* cosfourtheta + (cij_kl(4) + cij_kl(20))* sinphi* sinfourtheta) + &
+ sinphi* (cosfourtheta* ((cij_kl(15) + cij_kl(17))* cosphi + &
+ cij_kl(16)* sinphi) + (cij_kl(2) + cij_kl(7) - 2.d0* cij_kl(8) + cij_kl(21) + &
+ (cij_kl(2) - cij_kl(7) + cij_kl(21))* costwophi)* sinphi* sintwothetasq + &
+ (-cij_kl(13) + cij_kl(9)* sinphisq)* sinfourtheta) + &
+ cosphi* (8.d0* cij_kl(11)* costhetasq* sinphi*sinphisq* sinthetasq + &
+ (-cij_kl(14) + (cij_kl(10) + cij_kl(18))* sinphisq)* sinfourtheta)))
+
+ cij_kll(20) = ONE_EIGHTH* (2.d0* cosphi* costheta* (-3.d0* cij_kl(4) - cij_kl(9) + &
+ 4.d0* cij_kl(13) + cij_kl(20) + (cij_kl(4) + 3.d0* cij_kl(9) - &
+ 4.d0* cij_kl(13) + cij_kl(20))* costwotheta) + &
+ (cij_kl(4) - cij_kl(9) + cij_kl(20))* costhreephi* (costheta + &
+ 3.d0* costhreetheta) - &
+ 2.d0* costheta* (-cij_kl(5) - 3.d0* cij_kl(10) + 4.d0* cij_kl(14) + &
+ cij_kl(18) + (3.d0* cij_kl(5) + &
+ cij_kl(10) - 4.d0* cij_kl(14) + cij_kl(18))*costwotheta)* sinphi - &
+ (cij_kl(5) - cij_kl(10) - cij_kl(18))* &
+ (costheta + 3.d0* costhreetheta)* sinthreephi + 8.d0* (cij_kl(6) - &
+ cij_kl(11))* cosfourphi* costhetasq* sintheta - 8.d0* (cij_kl(1) - &
+ cij_kl(3) - cij_kl(7) + cij_kl(8) + &
+ (cij_kl(1) - cij_kl(3) - cij_kl(7) + cij_kl(8) + cij_kl(16) - &
+ cij_kl(19))* costwotheta)* sintwophi* sintheta - &
+ 8.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
+ cij_kl(21))* costhetasq* sinfourphi* sintheta + &
+ 2.d0* costwophi* ((cij_kl(6) + cij_kl(11) - 2.d0* cij_kl(15) + &
+ 2.d0* cij_kl(17))* sintheta + &
+ (cij_kl(6) + cij_kl(11) - 2.d0* (cij_kl(15) + cij_kl(17)))* sinthreetheta))
+
+ cij_kll(21) = ONE_FOURTH* (cij_kl(1) - cij_kl(2) + cij_kl(7) + cij_kl(16) + &
+ cij_kl(19) + cij_kl(21) - 2.d0* (cij_kl(1) - cij_kl(2) + cij_kl(7) - &
+ cij_kl(21))* cosfourphi* costhetasq + &
+ (cij_kl(1) - cij_kl(2) + cij_kl(7) - cij_kl(16) - cij_kl(19) + &
+ cij_kl(21))* costwotheta + &
+ 2.d0* (-cij_kl(6) + cij_kl(11))* costhetasq* sinfourphi - &
+ 2.d0* ((-cij_kl(16) + cij_kl(19))* costwophi + cij_kl(17)* sintwophi)* sinthetasq - &
+ ((cij_kl(5) - cij_kl(10) + cij_kl(18))* cosphi + (-cij_kl(5) + cij_kl(10) +&
+ cij_kl(18))* costhreephi + &
+ (-cij_kl(4) + cij_kl(9) + cij_kl(20))* sinphi - (cij_kl(4) - cij_kl(9) + &
+ cij_kl(20))* sinthreephi)* sintwotheta)
+
+ end subroutine rotate_kernels_dble
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/expand_compute_strain_product.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/expand_compute_strain_product.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/expand_compute_strain_product.f90 2013-07-13 14:08:06 UTC (rev 22584)
@@ -0,0 +1,25 @@
+
+ program expand_compute_strain_product
+
+!! DK DK July 2013: expand the calculations in compute_strain_product() in order to speed up the code
+
+ implicit none
+
+ integer :: i,j,p
+
+ ! Computing the 21 strain products without assuming eps(i)*b_eps(j) = eps(j)*b_eps(i)
+ p=1
+ do i=1,6
+ do j=i,6
+ print *,'prod(',p,')=eps',i,'*b_eps',j
+ if(j>i) then
+ print *,'prod(',p,')=prod(',p,')+eps',j,'*b_eps',i
+ if(j>3 .and. i<4) print *,'prod(',p,') = prod(',p,') * 2._CUSTOM_REAL'
+ endif
+ if(i>3) print *,'prod(',p,') = prod(',p,') * 4._CUSTOM_REAL'
+ p=p+1
+ enddo
+ enddo
+
+ end program expand_compute_strain_product
+
More information about the CIG-COMMITS
mailing list