[cig-commits] r21573 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Mar 19 14:31:52 PDT 2013
Author: dkomati1
Date: 2013-03-19 14:31:52 -0700 (Tue, 19 Mar 2013)
New Revision: 21573
Added:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_Dev.f90
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
Log:
prepared the framework for a future compute_forces_acoustic_Dev.f90 routine (i.e. adding USE_DEVILLE support for acoustic elements)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2013-03-19 21:15:27 UTC (rev 21572)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2013-03-19 21:31:52 UTC (rev 21573)
@@ -187,6 +187,7 @@
$O/compute_coupling_poroelastic_el.o \
$O/compute_forces_acoustic_calling_routine.o \
$O/compute_forces_acoustic_noDev.o \
+ $O/compute_forces_acoustic_Dev.o \
$O/compute_forces_viscoelastic_calling_routine.o \
$O/compute_forces_viscoelastic_Dev.o \
$O/compute_forces_viscoelastic_noDev.o \
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_Dev.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_Dev.f90 2013-03-19 21:31:52 UTC (rev 21573)
@@ -0,0 +1,184 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 2 . 1
+! ---------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA and CNRS / INRIA / University of Pau
+! (c) Princeton University / California Institute of Technology and CNRS / INRIA / University of Pau
+! July 2012
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! for acoustic solver
+
+ subroutine compute_forces_acoustic_Dev(iphase,NSPEC_AB,NGLOB_AB, &
+ potential_acoustic,potential_dot_dot_acoustic, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ rhostore,jacobian,ibool, &
+ num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
+ phase_ispec_inner_acoustic)
+
+! computes forces for acoustic elements
+!
+! note that pressure is defined as:
+! p = - Chi_dot_dot
+!
+ use specfem_par,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,TINYVAL_SNGL,ABSORB_USE_PML,ABSORBING_CONDITIONS,PML_CONDITIONS
+
+ implicit none
+
+ !include "constants.h"
+ integer :: NSPEC_AB,NGLOB_AB
+
+ ! acoustic potentials
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
+ potential_acoustic,potential_dot_dot_acoustic
+
+ ! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ rhostore,jacobian
+
+ ! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+ integer :: iphase
+ integer :: num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic
+ integer, dimension(num_phase_ispec_acoustic,2) :: phase_ispec_inner_acoustic
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: chi_elem
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: temp1,temp2,temp3
+ real(kind=CUSTOM_REAL) :: temp1l,temp2l,temp3l
+
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) :: dpotentialdxl,dpotentialdyl,dpotentialdzl
+ real(kind=CUSTOM_REAL) :: rho_invl
+
+ integer :: ispec,iglob,i,j,k,l,ispec_p,num_elements
+
+ if( iphase == 1 ) then
+ num_elements = nspec_outer_acoustic
+ else
+ num_elements = nspec_inner_acoustic
+ endif
+
+ ! loop over spectral elements
+ do ispec_p = 1,num_elements
+
+ ispec = phase_ispec_inner_acoustic(ispec_p,iphase)
+
+ ! gets values for element
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ chi_elem(i,j,k) = potential_acoustic(ibool(i,j,k,ispec))
+ enddo
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ ! derivative along x, y, z
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the loops because NGLLX == NGLLY == NGLLZ
+ temp1l = 0._CUSTOM_REAL
+ temp2l = 0._CUSTOM_REAL
+ temp3l = 0._CUSTOM_REAL
+
+ do l = 1,NGLLX
+ temp1l = temp1l + chi_elem(l,j,k)*hprime_xx(i,l)
+ temp2l = temp2l + chi_elem(i,l,k)*hprime_yy(j,l)
+ temp3l = temp3l + chi_elem(i,j,l)*hprime_zz(k,l)
+ enddo
+
+ ! get derivatives of potential with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
+
+ ! derivatives of potential
+ dpotentialdxl = xixl*temp1l + etaxl*temp2l + gammaxl*temp3l
+ dpotentialdyl = xiyl*temp1l + etayl*temp2l + gammayl*temp3l
+ dpotentialdzl = xizl*temp1l + etazl*temp2l + gammazl*temp3l
+
+ ! density (reciproc)
+ rho_invl = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec)
+
+ ! for acoustic medium
+ ! also add GLL integration weights
+ temp1(i,j,k) = rho_invl * wgllwgll_yz(j,k) * jacobianl* &
+ (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ temp2(i,j,k) = rho_invl * wgllwgll_xz(i,k) * jacobianl* &
+ (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ temp3(i,j,k) = rho_invl * wgllwgll_xy(i,j) * jacobianl* &
+ (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+ enddo
+ enddo
+ enddo
+
+ ! second double-loop over GLL to compute all the terms
+ do k = 1,NGLLZ
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+ ! along x,y,z direction
+ ! and assemble the contributions
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ temp1l = 0._CUSTOM_REAL
+ temp2l = 0._CUSTOM_REAL
+ temp3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ temp1l = temp1l + temp1(l,j,k) * hprimewgll_xx(l,i)
+ temp2l = temp2l + temp2(i,l,k) * hprimewgll_yy(l,j)
+ temp3l = temp3l + temp3(i,j,l) * hprimewgll_zz(l,k)
+ enddo
+
+ ! sum contributions from each element to the global values
+ iglob = ibool(i,j,k,ispec)
+
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - ( temp1l + temp2l + temp3l )
+
+ enddo
+ enddo
+ enddo
+
+ enddo ! end of loop over all spectral elements
+
+ end subroutine compute_forces_acoustic_Dev
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-03-19 21:15:27 UTC (rev 21572)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-03-19 21:31:52 UTC (rev 21573)
@@ -96,7 +96,19 @@
! acoustic pressure term
if(.NOT. GPU_MODE) then
! on CPU
- call compute_forces_acoustic_noDev( iphase, NSPEC_AB,NGLOB_AB, &
+ if(USE_DEVILLE_PRODUCTS) then
+ ! uses Deville (2002) optimizations
+ call compute_forces_acoustic_Dev(iphase,NSPEC_AB,NGLOB_AB, &
+ potential_acoustic,potential_dot_dot_acoustic, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ rhostore,jacobian,ibool, &
+ num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
+ phase_ispec_inner_acoustic)
+ else
+ call compute_forces_acoustic_noDev(iphase,NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz, &
@@ -106,11 +118,14 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic )
+ phase_ispec_inner_acoustic)
+ endif
! adjoint simulations
- if( SIMULATION_TYPE == 3 ) &
- call compute_forces_acoustic_noDev( iphase, NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ if( SIMULATION_TYPE == 3 ) then
+ if(USE_DEVILLE_PRODUCTS) then
+ ! uses Deville (2002) optimizations
+ call compute_forces_acoustic_Dev(iphase,NSPEC_ADJOINT,NGLOB_ADJOINT, &
b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz, &
@@ -120,7 +135,22 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic )
+ phase_ispec_inner_acoustic)
+ else
+ call compute_forces_acoustic_noDev(iphase,NSPEC_ADJOINT,NGLOB_ADJOINT, &
+ b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ rhostore,jacobian,ibool,deltat, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
+ phase_ispec_inner_acoustic)
+ endif
+ endif
+
else
! on GPU
! includes code for SIMULATION_TYPE==3
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-03-19 21:15:27 UTC (rev 21572)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-03-19 21:31:52 UTC (rev 21573)
@@ -26,7 +26,7 @@
! for acoustic solver
- subroutine compute_forces_acoustic_noDev( iphase, NSPEC_AB,NGLOB_AB, &
+ subroutine compute_forces_acoustic_noDev(iphase,NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz, &
@@ -36,7 +36,7 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic )
+ phase_ispec_inner_acoustic)
! computes forces for acoustic elements
!
@@ -74,12 +74,6 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
-! communication overlap
-! logical, dimension(NSPEC_AB) :: ispec_is_inner
-! logical :: phase_is_inner
-
-! logical, dimension(NSPEC_AB) :: ispec_is_acoustic
-
integer :: iphase
integer :: num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic
integer, dimension(num_phase_ispec_acoustic,2) :: phase_ispec_inner_acoustic
@@ -126,12 +120,8 @@
! loop over spectral elements
do ispec_p = 1,num_elements
- !if( (ispec_is_inner(ispec) .eqv. phase_is_inner) ) then
-
ispec = phase_ispec_inner_acoustic(ispec_p,iphase)
- !if( ispec_is_acoustic(ispec) ) then
-
! gets values for element
do k=1,NGLLZ
do j=1,NGLLY
@@ -141,9 +131,6 @@
enddo
enddo
- ! would check if anything to do, but might lower accuracy of computation
- !if( maxval( abs( chi_elem ) ) < TINYVAL_SNGL ) cycle
-
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -170,7 +157,6 @@
hp1 = hprime_xx(l,i)
iglob = ibool(l,j,k,ispec)
temp1l_new = temp1l_new + deltat*potential_dot_acoustic(iglob)*hp1
-
!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
@@ -279,9 +265,6 @@
enddo
enddo
-! endif ! end of test if acoustic element
-! endif ! ispec_is_inner
-
enddo ! end of loop over all spectral elements
! C-PML boundary
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-03-19 21:15:27 UTC (rev 21572)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-03-19 21:31:52 UTC (rev 21573)
@@ -2288,7 +2288,7 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
+
! compute stress sigma
sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
sigma_yx = mul*duxdyl_x + mul*duydxl_x
More information about the CIG-COMMITS
mailing list