[cig-commits] r11924 - seismo/3D/CPML/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu May 8 14:21:09 PDT 2008
Author: dkomati1
Date: 2008-05-08 14:21:09 -0700 (Thu, 08 May 2008)
New Revision: 11924
Modified:
seismo/3D/CPML/trunk/seismic_CPML_2D_aniso.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_iso_second.f90
seismo/3D/CPML/trunk/seismic_CPML_3D_iso_MPI_OpenMP.f90
seismo/3D/CPML/trunk/seismic_PML_Collino_2D_iso.f90
seismo/3D/CPML/trunk/seismic_PML_Collino_3D_iso_OpenMP.f90
Log:
fixed two bugs found by Roland Martin: K damping not correctly staggered for one variable of the 2D CPML codes, and upper bounds shifted by 1 for some PML layers.
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_aniso.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_aniso.f90 2008-05-07 13:37:59 UTC (rev 11923)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_aniso.f90 2008-05-08 21:21:09 UTC (rev 11924)
@@ -455,7 +455,7 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
yoriginbottom = thickness_PML_y
- yorigintop = NY*DELTAY - thickness_PML_y
+ yorigintop = (NY-1)*DELTAY - thickness_PML_y
do j = 1,NY
@@ -594,7 +594,7 @@
memory_dvx_dy(i,j) = b_y_half(j) * memory_dvx_dy(i,j) + a_y_half(j) * value_dvx_dy
value_dvy_dx = value_dvy_dx / K_x(i) + memory_dvy_dx(i,j)
- value_dvx_dy = value_dvx_dy / K_y(j) + memory_dvx_dy(i,j)
+ value_dvx_dy = value_dvx_dy / K_y_half(j) + memory_dvx_dy(i,j)
sigmaxy(i,j) = sigmaxy(i,j) + c33 * (value_dvy_dx + value_dvx_dy) * DELTAT
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_iso_second.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_iso_second.f90 2008-05-07 13:37:59 UTC (rev 11923)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_iso_second.f90 2008-05-08 21:21:09 UTC (rev 11924)
@@ -366,7 +366,7 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
yoriginbottom = thickness_PML_y
- yorigintop = NY*DELTAY - thickness_PML_y
+ yorigintop = (NY-1)*DELTAY - thickness_PML_y
do j = 1,NY
@@ -561,7 +561,7 @@
memory_dvx_dy(i,j) = b_y_half(j) * memory_dvx_dy(i,j) + a_y_half(j) * value_dvx_dy
value_dvy_dx = value_dvy_dx / K_x(i) + memory_dvy_dx(i,j)
- value_dvx_dy = value_dvx_dy / K_y(j) + memory_dvx_dy(i,j)
+ value_dvx_dy = value_dvx_dy / K_y_half(j) + memory_dvx_dy(i,j)
sigmaxy(i,j) = sigmaxy(i,j) + mu_half_y * (value_dvy_dx + value_dvx_dy) * DELTAT
@@ -661,16 +661,16 @@
! in order to interpolate density at the right location in the staggered grid cell
! but in a homogeneous medium we can safely ignore it
total_energy_kinetic(it) = 0.5d0 * sum( &
- rho(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1)*( &
- vx(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1)**2 + &
- vy(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1)**2))
+ rho(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML)*( &
+ vx(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML)**2 + &
+ vy(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML)**2))
! add potential energy, defined as 1/2 epsilon_ij sigma_ij
! in principle we should interpolate the medium parameters at the right location
! in the staggered grid cell but in a homogeneous medium we can safely ignore it
total_energy_potential(it) = ZERO
- do j = NPOINTS_PML, NY-NPOINTS_PML+1
- do i = NPOINTS_PML, NX-NPOINTS_PML+1
+ do j = NPOINTS_PML+1, NY-NPOINTS_PML
+ do i = NPOINTS_PML+1, NX-NPOINTS_PML
epsilon_xx = ((lambda(i,j) + 2.d0*mu(i,j)) * sigmaxx(i,j) - lambda(i,j) * &
sigmayy(i,j)) / (4.d0 * mu(i,j) * (lambda(i,j) + mu(i,j)))
epsilon_yy = ((lambda(i,j) + 2.d0*mu(i,j)) * sigmayy(i,j) - lambda(i,j) * &
Modified: seismo/3D/CPML/trunk/seismic_CPML_3D_iso_MPI_OpenMP.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_3D_iso_MPI_OpenMP.f90 2008-05-07 13:37:59 UTC (rev 11923)
+++ seismo/3D/CPML/trunk/seismic_CPML_3D_iso_MPI_OpenMP.f90 2008-05-08 21:21:09 UTC (rev 11924)
@@ -484,7 +484,7 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
yoriginbottom = thickness_PML_y
- yorigintop = NY*DELTAY - thickness_PML_y
+ yorigintop = (NY-1)*DELTAY - thickness_PML_y
do j = 1,NY
@@ -555,7 +555,7 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
zoriginbottom = thickness_PML_z
- zorigintop = NZ*DELTAZ - thickness_PML_z
+ zorigintop = (NZ-1)*DELTAZ - thickness_PML_z
do k = 1,NZ
@@ -1087,15 +1087,15 @@
kmin = 1
kmax = NZ_LOCAL
- if(rank == 0) kmin = NPOINTS_PML
- if(rank == nb_procs-1) kmax = NZ_LOCAL-NPOINTS_PML+1
+ if(rank == 0) kmin = NPOINTS_PML+1
+ if(rank == nb_procs-1) kmax = NZ_LOCAL-NPOINTS_PML
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k,epsilon_xx,epsilon_yy,epsilon_zz,epsilon_xy,epsilon_xz,epsilon_yz) &
!$OMP SHARED(kmin,kmax,vx,vy,vz,sigmaxx,sigmayy,sigmazz, &
!$OMP sigmaxy,sigmaxz,sigmayz) REDUCTION(+:total_energy_kinetic,total_energy_potential)
do k = kmin,kmax
- do j = NPOINTS_PML, NY-NPOINTS_PML+1
- do i = NPOINTS_PML, NX-NPOINTS_PML+1
+ do j = NPOINTS_PML+1, NY-NPOINTS_PML
+ do i = NPOINTS_PML+1, NX-NPOINTS_PML
! compute kinetic energy first, defined as 1/2 rho ||v||^2
! in principle we should use rho_half_x_half_y instead of rho for vy
Modified: seismo/3D/CPML/trunk/seismic_PML_Collino_2D_iso.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_PML_Collino_2D_iso.f90 2008-05-07 13:37:59 UTC (rev 11923)
+++ seismo/3D/CPML/trunk/seismic_PML_Collino_2D_iso.f90 2008-05-08 21:21:09 UTC (rev 11924)
@@ -217,11 +217,11 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
xoriginleft = delta
- xoriginright = NX*h - delta
+ xoriginright = (NX-1)*h - delta
do i=1,NX
- xval = h*dble(i)
+ xval = h*dble(i-1)
if(xval < xoriginleft) then
dx_over_two(i) = d0 * ((xoriginleft-xval)/delta)**2
@@ -243,11 +243,11 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
xoriginleft = delta
- xoriginright = NY*h - delta
+ xoriginright = (NY-1)*h - delta
do j=1,NY
- xval = h*dble(j)
+ xval = h*dble(j-1)
if(xval < xoriginleft) then
dy_over_two(j) = d0 * ((xoriginleft-xval)/delta)**2
@@ -502,17 +502,17 @@
! in order to interpolate density at the right location in the staggered grid cell
! but in a homogeneous medium we can safely ignore it
total_energy_kinetic(it) = 0.5d0 * sum(rho*( &
- (vx_1(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1) + &
- vx_2(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1))**2 + &
- (vy_1(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1) + &
- vy_2(NPOINTS_PML:NX-NPOINTS_PML+1,NPOINTS_PML:NY-NPOINTS_PML+1))**2))
+ (vx_1(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML) + &
+ vx_2(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML))**2 + &
+ (vy_1(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML) + &
+ vy_2(NPOINTS_PML+1:NX-NPOINTS_PML,NPOINTS_PML+1:NY-NPOINTS_PML))**2))
! add potential energy, defined as 1/2 epsilon_ij sigma_ij
! in principle we should interpolate the medium parameters at the right location
! in the staggered grid cell but in a homogeneous medium we can safely ignore it
total_energy_potential(it) = ZERO
- do j = NPOINTS_PML, NY-NPOINTS_PML+1
- do i = NPOINTS_PML, NX-NPOINTS_PML+1
+ do j = NPOINTS_PML+1, NY-NPOINTS_PML
+ do i = NPOINTS_PML+1, NX-NPOINTS_PML
! compute total field from split components
sigmaxx_total = sigmaxx_1(i,j) + sigmaxx_2(i,j)
Modified: seismo/3D/CPML/trunk/seismic_PML_Collino_3D_iso_OpenMP.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_PML_Collino_3D_iso_OpenMP.f90 2008-05-07 13:37:59 UTC (rev 11923)
+++ seismo/3D/CPML/trunk/seismic_PML_Collino_3D_iso_OpenMP.f90 2008-05-08 21:21:09 UTC (rev 11924)
@@ -218,11 +218,11 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
xoriginleft = delta
- xoriginright = NX*h - delta
+ xoriginright = (NX-1)*h - delta
do i=1,NX
- xval = h*dble(i)
+ xval = h*dble(i-1)
if(xval < xoriginleft) then
dx_over_two(i) = d0 * ((xoriginleft-xval)/delta)**2
@@ -244,11 +244,11 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
xoriginleft = delta
- xoriginright = NY*h - delta
+ xoriginright = (NY-1)*h - delta
do j=1,NY
- xval = h*dble(j)
+ xval = h*dble(j-1)
if(xval < xoriginleft) then
dy_over_two(j) = d0 * ((xoriginleft-xval)/delta)**2
@@ -270,11 +270,11 @@
! origin of the PML layer (position of right edge minus thickness, in meters)
xoriginleft = delta
- xoriginright = NZ*h - delta
+ xoriginright = (NZ-1)*h - delta
do k=1,NZ
- xval = h*dble(k)
+ xval = h*dble(k-1)
if(xval < xoriginleft) then
dz_over_two(k) = d0 * ((xoriginleft-xval)/delta)**2
@@ -765,9 +765,9 @@
!$OMP SHARED(vx_1,vx_2,vx_3,vy_1,vy_2,vy_3,vz_1,vz_2,vz_3,sigmaxx_1,sigmaxx_2, &
!$OMP sigmaxx_3,sigmayy_1,sigmayy_2,sigmayy_3,sigmazz_1,sigmazz_2,sigmazz_3, &
!$OMP sigmaxy_1,sigmaxy_2,sigmaxz_1,sigmaxz_3,sigmayz_2,sigmayz_3) REDUCTION(+:total_energy_kinetic,total_energy_potential)
- do k = NPOINTS_PML, NZ-NPOINTS_PML+1
- do j = NPOINTS_PML, NY-NPOINTS_PML+1
- do i = NPOINTS_PML, NX-NPOINTS_PML+1
+ do k = NPOINTS_PML+1, NZ-NPOINTS_PML
+ do j = NPOINTS_PML+1, NY-NPOINTS_PML
+ do i = NPOINTS_PML+1, NX-NPOINTS_PML
! compute kinetic energy first, defined as 1/2 rho ||v||^2
! in principle we should use rho_half_x_half_y instead of rho for vy
More information about the cig-commits
mailing list