[cig-commits] [commit] devel: fixes CPML when using Deville routine in fluids (da0f08c)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 27 02:53:42 PST 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/13ec3e35977ac71f4785cfd2493cb65d604d3e63...c03570160da2904c3100ec5e80b5068857972a33
>---------------------------------------------------------------
commit da0f08cdd2cf59c791c362ad4380d47217c5add0
Author: daniel peter <peterda at ethz.ch>
Date: Thu Nov 27 11:39:43 2014 +0100
fixes CPML when using Deville routine in fluids
>---------------------------------------------------------------
da0f08cdd2cf59c791c362ad4380d47217c5add0
src/specfem3D/compute_forces_acoustic_Dev.F90 | 223 +++++++++++++-----------
src/specfem3D/compute_forces_acoustic_noDev.f90 | 70 ++++----
src/specfem3D/pml_compute_memory_variables.f90 | 54 +++---
src/specfem3D/pml_output_VTKs.f90 | 5 +-
4 files changed, 182 insertions(+), 170 deletions(-)
diff --git a/src/specfem3D/compute_forces_acoustic_Dev.F90 b/src/specfem3D/compute_forces_acoustic_Dev.F90
index d9080c9..4a61127 100644
--- a/src/specfem3D/compute_forces_acoustic_Dev.F90
+++ b/src/specfem3D/compute_forces_acoustic_Dev.F90
@@ -43,44 +43,45 @@
!
use specfem_par,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ, &
m1,m2,NGLLCUBE,PML_CONDITIONS
+
use pml_par, only: is_CPML, spec_to_CPML, NSPEC_CPML, &
PML_dpotential_dxl,PML_dpotential_dyl,PML_dpotential_dzl,&
PML_dpotential_dxl_old,PML_dpotential_dyl_old,PML_dpotential_dzl_old,&
+ PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new,&
potential_dot_dot_acoustic_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl,&
- rmemory_dpotential_dzl,rmemory_potential_acoustic,potential_acoustic_old
+ rmemory_dpotential_dzl,rmemory_potential_acoustic, &
+ potential_acoustic_old,potential_acoustic_new
implicit none
- integer :: NSPEC_AB,NGLOB_AB
+ integer,intent(in) :: NSPEC_AB,NGLOB_AB
! acoustic potentials
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB),intent(inout) :: &
potential_acoustic,potential_dot_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) :: &
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: &
rhostore,jacobian
! array with derivatives of Lagrange polynomials and precalculated products
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprimewgll_xx,hprimewgll_xxT
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX),intent(in) :: hprime_xx,hprime_xxT
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX),intent(in) :: hprimewgll_xx,hprimewgll_xxT
- 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),intent(in) :: wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D
-! CPML adjoint
- logical :: backward_simulation
+ integer,intent(in) :: iphase
+ integer,intent(in) :: num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic
+ integer, dimension(num_phase_ispec_acoustic,2),intent(in) :: phase_ispec_inner_acoustic
- integer :: ispec_CPML
+ ! CPML adjoint
+ logical,intent(in) :: backward_simulation
+ ! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: tempx1,tempx2,tempx3
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: tempx1_old,tempx2_old,tempx3_old
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: newtempx1,newtempx2,newtempx3
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
@@ -91,31 +92,47 @@
! manually inline the calls to the Deville et al. (2002) routines
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: chi_elem
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: chi_elem_old
real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points
real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points
- real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points_old
- real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points_old
real(kind=CUSTOM_REAL), dimension(m1,m2) :: E1_m1_m2_5points
equivalence(chi_elem,B1_m1_m2_5points)
equivalence(tempx1,C1_m1_m2_5points)
- equivalence(chi_elem_old,B1_m1_m2_5points_old)
- equivalence(tempx1_old,C1_m1_m2_5points_old)
equivalence(newtempx1,E1_m1_m2_5points)
real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points
real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1_mxm_m2_m1_5points
- real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points_old
- real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1_mxm_m2_m1_5points_old
real(kind=CUSTOM_REAL), dimension(m2,m1) :: E1_mxm_m2_m1_5points
equivalence(chi_elem,A1_mxm_m2_m1_5points)
equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(newtempx3,E1_mxm_m2_m1_5points)
+
+ ! CPML
+ integer :: ispec_CPML
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: tempx1_old,tempx2_old,tempx3_old
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: tempx1_new,tempx2_new,tempx3_new
+
+ ! CPML Deville
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: chi_elem_old,chi_elem_new
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points_old,B1_m1_m2_5points_new
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points_old,C1_m1_m2_5points_new
+
+ equivalence(chi_elem_old,B1_m1_m2_5points_old)
+ equivalence(tempx1_old,C1_m1_m2_5points_old)
+
+ equivalence(chi_elem_new,B1_m1_m2_5points_new)
+ equivalence(tempx1_new,C1_m1_m2_5points_new)
+
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points_old,A1_mxm_m2_m1_5points_new
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1_mxm_m2_m1_5points_old,C1_mxm_m2_m1_5points_new
+
equivalence(chi_elem_old,A1_mxm_m2_m1_5points_old)
equivalence(tempx3_old,C1_mxm_m2_m1_5points_old)
- equivalence(newtempx3,E1_mxm_m2_m1_5points)
+
+ equivalence(chi_elem_new,A1_mxm_m2_m1_5points_new)
+ equivalence(tempx3_new,C1_mxm_m2_m1_5points_new)
#ifdef FORCE_VECTORIZATION
integer :: ijk
@@ -195,6 +212,7 @@
do j=1,NGLLY
do i=1,NGLLX
chi_elem_old(i,j,k) = potential_acoustic_old(ibool(i,j,k,ispec))
+ chi_elem_new(i,j,k) = potential_acoustic_new(ibool(i,j,k,ispec))
enddo
enddo
enddo
@@ -203,6 +221,7 @@
! thus use only for production runs with no bound checking
do ijk = 1,NGLLCUBE
chi_elem_old(ijk,1,1) = potential_acoustic_old(ibool(ijk,1,1,ispec))
+ chi_elem_new(ijk,1,1) = potential_acoustic_new(ibool(ijk,1,1,ispec))
enddo
#endif
@@ -216,6 +235,12 @@
hprime_xx(i,3)*B1_m1_m2_5points_old(3,j) + &
hprime_xx(i,4)*B1_m1_m2_5points_old(4,j) + &
hprime_xx(i,5)*B1_m1_m2_5points_old(5,j)
+
+ C1_m1_m2_5points_new(i,j) = hprime_xx(i,1)*B1_m1_m2_5points_new(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points_new(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points_new(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points_new(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points_new(5,j)
enddo
enddo
@@ -227,6 +252,12 @@
chi_elem_old(i,3,k)*hprime_xxT(3,j) + &
chi_elem_old(i,4,k)*hprime_xxT(4,j) + &
chi_elem_old(i,5,k)*hprime_xxT(5,j)
+
+ tempx2_new(i,j,k) = chi_elem_new(i,1,k)*hprime_xxT(1,j) + &
+ chi_elem_new(i,2,k)*hprime_xxT(2,j) + &
+ chi_elem_new(i,3,k)*hprime_xxT(3,j) + &
+ chi_elem_new(i,4,k)*hprime_xxT(4,j) + &
+ chi_elem_new(i,5,k)*hprime_xxT(5,j)
enddo
enddo
enddo
@@ -238,72 +269,22 @@
A1_mxm_m2_m1_5points_old(i,3)*hprime_xxT(3,j) + &
A1_mxm_m2_m1_5points_old(i,4)*hprime_xxT(4,j) + &
A1_mxm_m2_m1_5points_old(i,5)*hprime_xxT(5,j)
- enddo
- enddo
-
-#ifndef FORCE_VECTORIZATION
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ! 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
- PML_dpotential_dxl(i,j,k) = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
- PML_dpotential_dyl(i,j,k) = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
- PML_dpotential_dzl(i,j,k) = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
-
- PML_dpotential_dxl_old(i,j,k) = xixl*tempx1_old(i,j,k) + etaxl*tempx2_old(i,j,k) + gammaxl*tempx3_old(i,j,k)
- PML_dpotential_dyl_old(i,j,k) = xiyl*tempx1_old(i,j,k) + etayl*tempx2_old(i,j,k) + gammayl*tempx3_old(i,j,k)
- PML_dpotential_dzl_old(i,j,k) = xizl*tempx1_old(i,j,k) + etazl*tempx2_old(i,j,k) + gammazl*tempx3_old(i,j,k)
-
- enddo
+ C1_mxm_m2_m1_5points_new(i,j) = A1_mxm_m2_m1_5points_new(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points_new(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points_new(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points_new(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points_new(i,5)*hprime_xxT(5,j)
enddo
enddo
-#else
- do ijk = 1,NGLLCUBE
- ! get derivatives of potential with respect to x, y and z
- xixl = xix(ijk,1,1,ispec)
- xiyl = xiy(ijk,1,1,ispec)
- xizl = xiz(ijk,1,1,ispec)
- etaxl = etax(ijk,1,1,ispec)
- etayl = etay(ijk,1,1,ispec)
- etazl = etaz(ijk,1,1,ispec)
- gammaxl = gammax(ijk,1,1,ispec)
- gammayl = gammay(ijk,1,1,ispec)
- gammazl = gammaz(ijk,1,1,ispec)
- jacobianl = jacobian(ijk,1,1,ispec)
-
- ! derivatives of potential
- PML_dpotential_dxl(ijk,1,1) = xixl*tempx1(ijk,1,1) + etaxl*tempx2(ijk,1,1) + gammaxl*tempx3(ijk,1,1)
- PML_dpotential_dyl(ijk,1,1) = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
- PML_dpotential_dzl(ijk,1,1) = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
-
- PML_dpotential_dxl_old(ijk,1,1) = xixl*tempx1_old(ijk,1,1) + etaxl*tempx2_old(ijk,1,1) + gammaxl*tempx3_old(ijk,1,1)
- PML_dpotential_dyl_old(ijk,1,1) = xiyl*tempx1_old(ijk,1,1) + etayl*tempx2_old(ijk,1,1) + gammayl*tempx3_old(ijk,1,1)
- PML_dpotential_dzl_old(ijk,1,1) = xizl*tempx1_old(ijk,1,1) + etazl*tempx2_old(ijk,1,1) + gammazl*tempx3_old(ijk,1,1)
-
- enddo
-#endif
- endif
- endif
+ endif ! is_CPML
+ endif ! PML_CONDITIONS
#ifndef FORCE_VECTORIZATION
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
-
! get derivatives of potential with respect to x, y and z
xixl = xix(i,j,k,ispec)
xiyl = xiy(i,j,k,ispec)
@@ -321,17 +302,33 @@
dpotentialdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
dpotentialdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+ ! stores derivatives of ux, uy and uz with respect to x, y and z
+ if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array is_CPML() is unallocated when PML_CONDITIONS is false
+ if (is_CPML(ispec)) then
+ PML_dpotential_dxl(i,j,k) = dpotentialdxl
+ PML_dpotential_dyl(i,j,k) = dpotentialdyl
+ PML_dpotential_dzl(i,j,k) = dpotentialdzl
+
+ PML_dpotential_dxl_old(i,j,k) = xixl*tempx1_old(i,j,k) + etaxl*tempx2_old(i,j,k) + gammaxl*tempx3_old(i,j,k)
+ PML_dpotential_dyl_old(i,j,k) = xiyl*tempx1_old(i,j,k) + etayl*tempx2_old(i,j,k) + gammayl*tempx3_old(i,j,k)
+ PML_dpotential_dzl_old(i,j,k) = xizl*tempx1_old(i,j,k) + etazl*tempx2_old(i,j,k) + gammazl*tempx3_old(i,j,k)
+
+ PML_dpotential_dxl_new(i,j,k) = xixl*tempx1_new(i,j,k) + etaxl*tempx2_new(i,j,k) + gammaxl*tempx3_new(i,j,k)
+ PML_dpotential_dyl_new(i,j,k) = xiyl*tempx1_new(i,j,k) + etayl*tempx2_new(i,j,k) + gammayl*tempx3_new(i,j,k)
+ PML_dpotential_dzl_new(i,j,k) = xizl*tempx1_new(i,j,k) + etazl*tempx2_new(i,j,k) + gammazl*tempx3_new(i,j,k)
+ endif
+ endif
+
! density (reciproc)
rho_invl = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec)
! for acoustic medium
! also add GLL integration weights
- tempx1(i,j,k) = rho_invl * jacobianl* &
- (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- tempx2(i,j,k) = rho_invl * jacobianl* &
- (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- tempx3(i,j,k) = rho_invl * jacobianl* &
- (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+ tempx1(i,j,k) = rho_invl * jacobianl * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ tempx2(i,j,k) = rho_invl * jacobianl * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ tempx3(i,j,k) = rho_invl * jacobianl * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
enddo
enddo
enddo
@@ -354,17 +351,33 @@
dpotentialdyl = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
dpotentialdzl = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
+ ! stores derivatives of ux, uy and uz with respect to x, y and z
+ if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array is_CPML() is unallocated when PML_CONDITIONS is false
+ if (is_CPML(ispec)) then
+ PML_dpotential_dxl(ijk,1,1) = dpotentialdxl
+ PML_dpotential_dyl(ijk,1,1) = dpotentialdyl
+ PML_dpotential_dzl(ijk,1,1) = dpotentialdzl
+
+ PML_dpotential_dxl_old(ijk,1,1) = xixl*tempx1_old(ijk,1,1) + etaxl*tempx2_old(ijk,1,1) + gammaxl*tempx3_old(ijk,1,1)
+ PML_dpotential_dyl_old(ijk,1,1) = xiyl*tempx1_old(ijk,1,1) + etayl*tempx2_old(ijk,1,1) + gammayl*tempx3_old(ijk,1,1)
+ PML_dpotential_dzl_old(ijk,1,1) = xizl*tempx1_old(ijk,1,1) + etazl*tempx2_old(ijk,1,1) + gammazl*tempx3_old(ijk,1,1)
+
+ PML_dpotential_dxl_new(ijk,1,1) = xixl*tempx1_new(ijk,1,1) + etaxl*tempx2_new(ijk,1,1) + gammaxl*tempx3_new(ijk,1,1)
+ PML_dpotential_dyl_new(ijk,1,1) = xiyl*tempx1_new(ijk,1,1) + etayl*tempx2_new(ijk,1,1) + gammayl*tempx3_new(ijk,1,1)
+ PML_dpotential_dzl_new(ijk,1,1) = xizl*tempx1_new(ijk,1,1) + etazl*tempx2_new(ijk,1,1) + gammazl*tempx3_new(ijk,1,1)
+ endif
+ endif
+
! density (reciproc)
rho_invl = 1.0_CUSTOM_REAL / rhostore(ijk,1,1,ispec)
! for acoustic medium
! also add GLL integration weights
- tempx1(ijk,1,1) = rho_invl * jacobianl* &
- (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- tempx2(ijk,1,1) = rho_invl * jacobianl* &
- (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- tempx3(ijk,1,1) = rho_invl * jacobianl* &
- (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+ tempx1(ijk,1,1) = rho_invl * jacobianl * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ tempx2(ijk,1,1) = rho_invl * jacobianl * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ tempx3(ijk,1,1) = rho_invl * jacobianl * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
enddo
#endif
@@ -372,8 +385,8 @@
! do not merge this second line with the first using an ".and." statement
! because array is_CPML() is unallocated when PML_CONDITIONS is false
if (is_CPML(ispec)) then
- tempx1 = 0._CUSTOM_REAL; tempx2 = 0._CUSTOM_REAL; tempx3 = 0._CUSTOM_REAL
ispec_CPML = spec_to_CPML(ispec)
+
! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
call pml_compute_memory_variables_acoustic(ispec,ispec_CPML,tempx1,tempx2,tempx3,&
rmemory_dpotential_dxl,rmemory_dpotential_dyl,rmemory_dpotential_dzl)
@@ -382,7 +395,7 @@
call pml_compute_accel_contribution_acoustic(ispec,ispec_CPML,potential_acoustic,&
potential_dot_acoustic,rmemory_potential_acoustic)
endif
- endif
+ endif ! PML_CONDITIONS
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
@@ -427,8 +440,10 @@
! 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) - (wgllwgll_yz_3D(i,j,k)*newtempx1(i,j,k) &
- + wgllwgll_xz_3D(i,j,k)*newtempx2(i,j,k) + wgllwgll_xy_3D(i,j,k)*newtempx3(i,j,k))
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - ( wgllwgll_yz_3D(i,j,k)*newtempx1(i,j,k) &
+ + wgllwgll_xz_3D(i,j,k)*newtempx2(i,j,k) &
+ + wgllwgll_xy_3D(i,j,k)*newtempx3(i,j,k))
enddo
enddo
@@ -441,8 +456,10 @@
do ijk = 1,NGLLCUBE
! sum contributions from each element to the global values
iglob = ibool(ijk,1,1,ispec)
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - (wgllwgll_yz_3D(ijk,1,1)*newtempx1(ijk,1,1) &
- + wgllwgll_xz_3D(ijk,1,1)*newtempx2(ijk,1,1) + wgllwgll_xy_3D(ijk,1,1)*newtempx3(ijk,1,1))
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - ( wgllwgll_yz_3D(ijk,1,1)*newtempx1(ijk,1,1) &
+ + wgllwgll_xz_3D(ijk,1,1)*newtempx2(ijk,1,1) &
+ + wgllwgll_xy_3D(ijk,1,1)*newtempx3(ijk,1,1))
enddo
#endif
@@ -456,8 +473,7 @@
do j = 1,NGLLZ
do i = 1,NGLLX
iglob = ibool(i,j,k,ispec)
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
- potential_dot_dot_acoustic_CPML(i,j,k)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_dot_acoustic_CPML(i,j,k)
enddo
enddo
enddo
@@ -469,8 +485,7 @@
enddo
#endif
endif
- endif
-
+ endif ! PML_CONDITIONS
enddo ! end of loop over all spectral elements
diff --git a/src/specfem3D/compute_forces_acoustic_noDev.f90 b/src/specfem3D/compute_forces_acoustic_noDev.f90
index 591c2a4..751a123 100644
--- a/src/specfem3D/compute_forces_acoustic_noDev.f90
+++ b/src/specfem3D/compute_forces_acoustic_noDev.f90
@@ -43,50 +43,54 @@
! p = - Chi_dot_dot
!
use specfem_par,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,PML_CONDITIONS
+
use pml_par, only: is_CPML, spec_to_CPML, NSPEC_CPML, &
PML_dpotential_dxl,PML_dpotential_dyl,PML_dpotential_dzl,&
PML_dpotential_dxl_old,PML_dpotential_dyl_old,PML_dpotential_dzl_old,&
PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new,&
potential_dot_dot_acoustic_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl,&
- rmemory_dpotential_dzl,rmemory_potential_acoustic,potential_acoustic_old,potential_acoustic_new
+ rmemory_dpotential_dzl,rmemory_potential_acoustic, &
+ potential_acoustic_old,potential_acoustic_new
implicit none
- integer :: NSPEC_AB,NGLOB_AB
+ integer,intent(in) :: NSPEC_AB,NGLOB_AB
! acoustic potentials
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB),intent(inout) :: &
potential_acoustic,potential_dot_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) :: &
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: &
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
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX),intent(in) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY),intent(in) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ),intent(in) :: hprime_zz,hprimewgll_zz
- 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),intent(in) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ),intent(in) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ),intent(in) :: wgllwgll_yz
-! CPML adjoint
- logical :: backward_simulation
+ integer,intent(in) :: iphase
+ integer,intent(in) :: num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic
+ integer, dimension(num_phase_ispec_acoustic,2),intent(in) :: phase_ispec_inner_acoustic
-! local variables
- real(kind=CUSTOM_REAL) :: hp1,hp2,hp3
+ ! CPML adjoint
+ logical,intent(in) :: backward_simulation
+ ! local variables
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) :: hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) :: fac1,fac2,fac3
+
+ ! CPML
real(kind=CUSTOM_REAL) :: temp1l_old,temp2l_old,temp3l_old
real(kind=CUSTOM_REAL) :: temp1l_new,temp2l_new,temp3l_new
@@ -203,7 +207,6 @@
PML_dpotential_dxl_new(i,j,k) = xixl*temp1l_new + etaxl*temp2l_new + gammaxl*temp3l_new
PML_dpotential_dyl_new(i,j,k) = xiyl*temp1l_new + etayl*temp2l_new + gammayl*temp3l_new
PML_dpotential_dzl_new(i,j,k) = xizl*temp1l_new + etazl*temp2l_new + gammazl*temp3l_new
-
endif
endif
@@ -211,13 +214,9 @@
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)
+ temp1(i,j,k) = rho_invl * jacobianl * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ temp2(i,j,k) = rho_invl * jacobianl * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ temp3(i,j,k) = rho_invl * jacobianl * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
enddo
enddo
enddo
@@ -242,7 +241,6 @@
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
-
! along x,y,z direction
! and assemble the contributions
!!! can merge these loops because NGLLX = NGLLY = NGLLZ
@@ -256,11 +254,14 @@
temp3l = temp3l + temp3(i,j,l) * hprimewgll_zz(l,k)
enddo
+ ! also add GLL integration weights
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
! 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 )
-
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - ( fac1*temp1l + fac2*temp2l + fac3*temp3l )
enddo
enddo
enddo
@@ -269,14 +270,11 @@
! do not merge this second line with the first using an ".and." statement
! because array is_CPML() is unallocated when PML_CONDITIONS is false
if (is_CPML(ispec)) then
-
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
iglob = ibool(i,j,k,ispec)
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
- potential_dot_dot_acoustic_CPML(i,j,k)
-
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_dot_acoustic_CPML(i,j,k)
enddo
enddo
enddo
diff --git a/src/specfem3D/pml_compute_memory_variables.f90 b/src/specfem3D/pml_compute_memory_variables.f90
index 3261b9e..18ba613 100644
--- a/src/specfem3D/pml_compute_memory_variables.f90
+++ b/src/specfem3D/pml_compute_memory_variables.f90
@@ -364,26 +364,28 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
use specfem_par, only: wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian,&
it,deltat,rhostore
+
use pml_par, only: NSPEC_CPML,CPML_regions,k_store_x,k_store_y,k_store_z,&
d_store_x,d_store_y,d_store_z,&
alpha_store_x,alpha_store_y,alpha_store_z,&
PML_dpotential_dxl,PML_dpotential_dyl,PML_dpotential_dzl,&
PML_dpotential_dxl_old,PML_dpotential_dyl_old,PML_dpotential_dzl_old,&
PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new
- use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,USE_DEVILLE_PRODUCTS, &
+
+ use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,FOUR_THIRDS, &
CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ
implicit none
integer, intent(in) :: ispec,ispec_CPML
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ), intent(out) :: temp1,temp2,temp3
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3) :: &
- rmemory_dpotential_dxl, rmemory_dpotential_dyl, rmemory_dpotential_dzl
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3),intent(inout) :: &
+ rmemory_dpotential_dxl, rmemory_dpotential_dyl, rmemory_dpotential_dzl
! local parameters
integer :: i,j,k
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
- real(kind=CUSTOM_REAL) :: rho_invl_jacob,rhoin_jacob_jk,rhoin_jacob_ik,rhoin_jacob_ij
+ real(kind=CUSTOM_REAL) :: rho_invl_jacob
real(kind=CUSTOM_REAL) :: dpotentialdxl,dpotentialdyl,dpotentialdzl
real(kind=CUSTOM_REAL) :: time_nplus1,time_n
real(kind=CUSTOM_REAL) :: A6,A7,A8,A9 ! L231
@@ -401,27 +403,6 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
- 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)
- rho_invl_jacob = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec) * jacobianl
- if (USE_DEVILLE_PRODUCTS) then
- rhoin_jacob_jk = rho_invl_jacob
- rhoin_jacob_ik = rho_invl_jacob
- rhoin_jacob_ij = rho_invl_jacob
- else
- rhoin_jacob_jk = rho_invl_jacob * wgllwgll_yz(j,k)
- rhoin_jacob_ik = rho_invl_jacob * wgllwgll_xz(i,k)
- rhoin_jacob_ij = rho_invl_jacob * wgllwgll_xy(i,j)
- endif
-
!---------------------- A6, A7, A8, A9 --------------------------
kappa_x = k_store_x(i,j,k,ispec_CPML)
kappa_y = k_store_y(i,j,k,ispec_CPML)
@@ -432,13 +413,13 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
alpha_x = alpha_store_x(i,j,k,ispec_CPML)
alpha_y = alpha_store_y(i,j,k,ispec_CPML)
alpha_z = alpha_store_z(i,j,k,ispec_CPML)
+
call lijk_parameter_computation(time_nplus1,deltat,&
kappa_z,d_z,alpha_z,kappa_y,d_y,alpha_y,kappa_x,d_x,alpha_x,&
CPML_region_local,231,A6,A7,A8,A9,&
coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2,&
coef0_3,coef1_3,coef2_3,singularity_type_2,singularity_type_3)
-
rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + &
coef1_1 * PML_dpotential_dxl_new(i,j,k) + coef2_1 * PML_dpotential_dxl_old(i,j,k)
@@ -458,6 +439,7 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
alpha_x = alpha_store_x(i,j,k,ispec_CPML)
alpha_y = alpha_store_y(i,j,k,ispec_CPML)
alpha_z = alpha_store_z(i,j,k,ispec_CPML)
+
call lijk_parameter_computation(time_nplus1,deltat,&
kappa_x,d_x,alpha_x,kappa_z,d_z,alpha_z,kappa_y,d_y,alpha_y,&
CPML_region_local,132,A10,A11,A12,A13,&
@@ -483,6 +465,7 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
alpha_x = alpha_store_x(i,j,k,ispec_CPML)
alpha_y = alpha_store_y(i,j,k,ispec_CPML)
alpha_z = alpha_store_z(i,j,k,ispec_CPML)
+
call lijk_parameter_computation(time_nplus1,deltat,&
kappa_x,d_x,alpha_x,kappa_y,d_y,alpha_y,kappa_z,d_z,alpha_z,&
CPML_region_local,123,A14,A15,A16,A17,&
@@ -510,10 +493,23 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te
A15 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + &
A16 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) + &
A17 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,3)
- temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+ 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)
+ rho_invl_jacob = jacobianl / rhostore(i,j,k,ispec)
+
+ temp1(i,j,k) = rho_invl_jacob * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ temp2(i,j,k) = rho_invl_jacob * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ temp3(i,j,k) = rho_invl_jacob * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
enddo
enddo
enddo
diff --git a/src/specfem3D/pml_output_VTKs.f90 b/src/specfem3D/pml_output_VTKs.f90
index caf7f6d..2c55cec 100644
--- a/src/specfem3D/pml_output_VTKs.f90
+++ b/src/specfem3D/pml_output_VTKs.f90
@@ -43,7 +43,10 @@ subroutine pml_output_VTKs()
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable:: temp_d_store_x,temp_d_store_y,temp_d_store_z
character(len=MAX_STRING_LEN) :: vtkfilename
- if (myrank == 0) write(IMAIN,*) 'Writing informations about C-PML elements in VTK-file format'
+ if (myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Writing informations about C-PML elements in VTK-file format'
+ endif
! C-PML regions
allocate(temp_CPML_regions(NSPEC_AB),stat=ier)
More information about the CIG-COMMITS
mailing list