[cig-commits] r18484 - seismo/3D/CPML/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat May 28 17:10:34 PDT 2011
Author: dkomati1
Date: 2011-05-28 17:10:34 -0700 (Sat, 28 May 2011)
New Revision: 18484
Modified:
seismo/3D/CPML/trunk/seismic_ADEPML_2D_RK4_eighth_order.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_second_order.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90
seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90
seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90
Log:
changed alpha_prime to alpha in variable names
Modified: seismo/3D/CPML/trunk/seismic_ADEPML_2D_RK4_eighth_order.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_ADEPML_2D_RK4_eighth_order.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_ADEPML_2D_RK4_eighth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -327,16 +327,16 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(-4:NX+4) :: d_x_1,K_x_1,alpha_prime_x_1,g_x_1,ksi_x
- double precision, dimension(-4:NX+4) :: d_x_half_1,K_x_half_1,alpha_prime_x_half_1,g_x_half_1,ksi_x_half
- double precision, dimension(-4:NY+4) :: d_y_1,K_y_1,alpha_prime_y_1,g_y_1,ksi_y
- double precision, dimension(-4:NY+4) :: d_y_half_1,K_y_half_1,alpha_prime_y_half_1,g_y_half_1,ksi_y_half
+ double precision, dimension(-4:NX+4) :: d_x_1,K_x_1,alpha_x_1,g_x_1,ksi_x
+ double precision, dimension(-4:NX+4) :: d_x_half_1,K_x_half_1,alpha_x_half_1,g_x_half_1,ksi_x_half
+ double precision, dimension(-4:NY+4) :: d_y_1,K_y_1,alpha_y_1,g_y_1,ksi_y
+ double precision, dimension(-4:NY+4) :: d_y_half_1,K_y_half_1,alpha_y_half_1,g_y_half_1,ksi_y_half
! 1D arrays for the damping profiles
- double precision, dimension(-4:NX+4) :: d_x_2,K_x_2,alpha_prime_x_2,g_x_2
- double precision, dimension(-4:NX+4) :: d_x_half_2,K_x_half_2,alpha_prime_x_half_2,g_x_half_2
- double precision, dimension(-4:NY+4) :: d_y_2,K_y_2,alpha_prime_y_2,g_y_2
- double precision, dimension(-4:NY+4) :: d_y_half_2,K_y_half_2,alpha_prime_y_half_2,g_y_half_2
+ double precision, dimension(-4:NX+4) :: d_x_2,K_x_2,alpha_x_2,g_x_2
+ double precision, dimension(-4:NX+4) :: d_x_half_2,K_x_half_2,alpha_x_half_2,g_x_half_2
+ double precision, dimension(-4:NY+4) :: d_y_2,K_y_2,alpha_y_2,g_y_2
+ double precision, dimension(-4:NY+4) :: d_y_half_2,K_y_half_2,alpha_y_half_2,g_y_half_2
! coefficients that allow to reset the memory variables at each RK4 substep depend on the substepping and are then of dimension 4,
! 1D arrays for the damping profiles
@@ -440,8 +440,8 @@
d_x_half_1(:) = ZERO
K_x_1(:) = 1.d0
K_x_half_1(:) = 1.d0
- alpha_prime_x_1(:) = ZERO
- alpha_prime_x_half_1(:) = ZERO
+ alpha_x_1(:) = ZERO
+ alpha_x_half_1(:) = ZERO
a_x_1(:,:) = ZERO
a_x_half_1(:,:) = ZERO
g_x_1(:) = 5.d-1
@@ -453,8 +453,8 @@
d_y_half_1(:) = ZERO
K_y_1(:) = 1.d0
K_y_half_1(:) = 1.d0
- alpha_prime_y_1(:) = ZERO
- alpha_prime_y_half_1(:) = ZERO
+ alpha_y_1(:) = ZERO
+ alpha_y_half_1(:) = ZERO
a_y_1(:,:) = ZERO
a_y_half_1(:,:) = ZERO
g_y_1(:) = 1.d0
@@ -464,8 +464,8 @@
d_x_half_2(:) = ZERO
K_x_2(:) = 1.d0
K_x_half_2(:) = 1.d0
- alpha_prime_x_2(:) = ZERO
- alpha_prime_x_half_2(:) = ZERO
+ alpha_x_2(:) = ZERO
+ alpha_x_half_2(:) = ZERO
a_x_2(:,:) = ZERO
a_x_half_2(:,:) = ZERO
g_x_2(:) = 1.d0
@@ -475,8 +475,8 @@
d_y_half_2(:) = ZERO
K_y_2(:) = 1.d0
K_y_half_2(:) = 1.d0
- alpha_prime_y_2(:) = ZERO
- alpha_prime_y_half_2(:) = ZERO
+ alpha_y_2(:) = ZERO
+ alpha_y_half_2(:) = ZERO
a_y_2(:,:) = ZERO
a_y_half_2(:,:) = ZERO
g_y_2(:) = 1.d0
@@ -503,7 +503,7 @@
d_x_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -513,7 +513,7 @@
d_x_half_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
@@ -528,7 +528,7 @@
d_x_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -538,7 +538,7 @@
d_x_half_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
@@ -548,29 +548,29 @@
d_x_half_2(i) = 0.d0
! just in case, for -5 at the end
- if(alpha_prime_x_1(i) < ZERO) alpha_prime_x_1(i) = ZERO
- if(alpha_prime_x_half_1(i) < ZERO) alpha_prime_x_half_1(i) = ZERO
+ if(alpha_x_1(i) < ZERO) alpha_x_1(i) = ZERO
+ if(alpha_x_half_1(i) < ZERO) alpha_x_half_1(i) = ZERO
! just in case, for -5 at the end
- if(alpha_prime_x_2(i) < ZERO) alpha_prime_x_2(i) = ZERO
- if(alpha_prime_x_half_2(i) < ZERO) alpha_prime_x_half_2(i) = ZERO
+ if(alpha_x_2(i) < ZERO) alpha_x_2(i) = ZERO
+ if(alpha_x_half_2(i) < ZERO) alpha_x_half_2(i) = ZERO
! CPML damping parameters for the 4 sub time steps of RK4 algorithm
do inc=1,4
- b_x_1(inc,i) = (1.-epsn*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_prime_x_1(i)))/&
- (1.+epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_prime_x_1(i)))
+ b_x_1(inc,i) = (1.-epsn*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_x_1(i)))/&
+ (1.+epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_x_1(i)))
b_x_half_1(inc,i) = (1.-epsn*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i) &
- + alpha_prime_x_half_1(i)))/(1. +epsn1*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i) &
- + alpha_prime_x_half_1(i)))
+ + alpha_x_half_1(i)))/(1. +epsn1*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i) &
+ + alpha_x_half_1(i)))
! this to avoid division by zero outside the PML
if(abs(d_x_1(i)) > 1.d-6) a_x_1(inc,i) = - DELTAT*rk41(inc)*d_x_1(i) / (K_x_1(i)* K_x_1(i))/&
- (1. +epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_prime_x_1(i)))
+ (1. +epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_x_1(i)))
if(abs(d_x_half_1(i)) > 1.d-6) a_x_half_1(inc,i) =-DELTAT*rk41(inc)*d_x_half_1(i)/&
(K_x_half_1(i)*K_x_half_1(i) )/&
(1. +epsn1*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i)&
- + alpha_prime_x_half_1(i)))
+ + alpha_x_half_1(i)))
enddo
@@ -597,7 +597,7 @@
d_y_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -607,7 +607,7 @@
d_y_half_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
@@ -622,7 +622,7 @@
d_y_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -632,30 +632,30 @@
d_y_half_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_y_1(j) < ZERO) alpha_prime_y_1(j) = ZERO
- if(alpha_prime_y_half_1(j) < ZERO) alpha_prime_y_half_1(j) = ZERO
+ if(alpha_y_1(j) < ZERO) alpha_y_1(j) = ZERO
+ if(alpha_y_half_1(j) < ZERO) alpha_y_half_1(j) = ZERO
! CPML damping parameters for the 4 sub time steps of RK4 algorithm
do inc=1,4
- b_y_1(inc,j) = (1.-epsn*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_prime_y_1(j)))/&
- (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_prime_y_1(j)))
+ b_y_1(inc,j) = (1.-epsn*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_y_1(j)))/&
+ (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_y_1(j)))
b_y_half_1(inc,j) = (1.-epsn*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j) + &
- alpha_prime_y_half_1(j)))/(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j)&
- + alpha_prime_y_half_1(j)))
+ alpha_y_half_1(j)))/(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j)&
+ + alpha_y_half_1(j)))
! this to avoid division by zero outside the PML
if(abs(d_y_1(j)) > 1.d-6) a_y_1(inc,j) = - DELTAT*rk41(inc)*d_y_1(j)&
/ (K_y_1(j)* K_y_1(j))/&
- (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_prime_y_1(j)))
+ (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_y_1(j)))
if(abs(d_y_half_1(j)) > 1.d-6) a_y_half_1(inc,j) = -DELTAT*rk41(inc)*d_y_half_1(j) /&
(K_y_half_1(j) * K_y_half_1(j) )/&
-(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j) + alpha_prime_y_half_1(j)))
+(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j) + alpha_y_half_1(j)))
enddo
enddo
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -279,7 +279,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -304,8 +304,8 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
double precision :: thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision :: Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -397,8 +397,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -406,8 +406,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -432,7 +432,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -442,7 +442,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -457,7 +457,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -467,22 +467,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -507,7 +507,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -517,7 +517,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -532,7 +532,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -542,18 +542,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -221,7 +221,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -246,8 +246,8 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
double precision :: thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision :: Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -310,8 +310,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -319,8 +319,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -345,7 +345,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -355,7 +355,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -370,7 +370,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -380,22 +380,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -420,7 +420,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -430,7 +430,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -445,7 +445,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -455,18 +455,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_second_order.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_second_order.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_second_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -220,7 +220,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -245,8 +245,8 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
double precision :: thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision :: Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -309,8 +309,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -318,8 +318,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -344,7 +344,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -354,7 +354,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -369,7 +369,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -379,22 +379,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -419,7 +419,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -429,7 +429,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -444,7 +444,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -454,18 +454,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -292,7 +292,7 @@
! double precision, parameter :: K_MAX_PML = 7.d0 ! from Gedney page 8.11
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
! double precision, parameter :: ALPHA_MAX_PML = 0.05d0 ! from Gedney page 8.22
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! 2D arrays for the memory variables
double precision, dimension(0:NX+1,0:NY+1) :: gamma11,gamma22
@@ -306,8 +306,8 @@
memory_dy_sigmaxx,memory_dy_sigmayy,memory_dy_sigmaxy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half_x,K_x_half_x,alpha_prime_x_half_x,a_x_half_x,b_x_half_x
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half_y,K_y_half_y,alpha_prime_y_half_y,a_y_half_y,b_y_half_y
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half_x,K_x_half_x,alpha_x_half_x,a_x_half_x,b_x_half_x
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half_y,K_y_half_y,alpha_y_half_y,a_y_half_y,b_y_half_y
double precision thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -421,11 +421,11 @@
K_y(:) = 1.d0
K_y_half_y(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half_x(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half_x(:) = ZERO
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half_y(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half_y(:) = ZERO
a_x(:) = ZERO
a_x_half_x(:) = ZERO
@@ -452,7 +452,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -462,7 +462,7 @@
d_x_half_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -477,7 +477,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -487,23 +487,23 @@
d_x_half_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half_x(i) < ZERO) alpha_prime_x_half_x(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half_x(i) < ZERO) alpha_x_half_x(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half_x(i) = exp(- (d_x_half_x(i) / K_x_half_x(i) + alpha_prime_x_half_x(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half_x(i) = exp(- (d_x_half_x(i) / K_x_half_x(i) + alpha_x_half_x(i)) * DELTAT)
! this to avoid division by zero outside the PML
if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) /&
- (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half_x(i)) > 1.d-6) a_x_half_x(i) = d_x_half_x(i)&
- * (b_x_half_x(i) - 1.d0) / (K_x_half_x(i) * (d_x_half_x(i) + K_x_half_x(i) * alpha_prime_x_half_x(i)))
+ * (b_x_half_x(i) - 1.d0) / (K_x_half_x(i) * (d_x_half_x(i) + K_x_half_x(i) * alpha_x_half_x(i)))
enddo
@@ -528,7 +528,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -538,7 +538,7 @@
d_y_half_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -553,7 +553,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -563,23 +563,23 @@
d_y_half_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
-! if(alpha_prime_y(j) < ZERO) alpha_prime_y(j) = ZERO
-! if(alpha_prime_y_half_y(j) < ZERO) alpha_prime_y_half_y(j) = ZERO
+! if(alpha_y(j) < ZERO) alpha_y(j) = ZERO
+! if(alpha_y_half_y(j) < ZERO) alpha_y_half_y(j) = ZERO
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half_y(j) = exp(- (d_y_half_y(j) / K_y_half_y(j) + alpha_prime_y_half_y(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half_y(j) = exp(- (d_y_half_y(j) / K_y_half_y(j) + alpha_y_half_y(j)) * DELTAT)
! this to avoid division by zero outside the PML
if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) &
- / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half_y(j)) > 1.d-6) a_y_half_y(j) = d_y_half_y(j)&
- * (b_y_half_y(j) - 1.d0) / (K_y_half_y(j) * (d_y_half_y(j) + K_y_half_y(j) * alpha_prime_y_half_y(j)))
+ * (b_y_half_y(j) - 1.d0) / (K_y_half_y(j) * (d_y_half_y(j) + K_y_half_y(j) * alpha_y_half_y(j)))
enddo
Modified: seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -208,7 +208,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -253,9 +253,9 @@
value_dsigmayz_dz
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
- double precision, dimension(NZ) :: d_z,K_z,alpha_prime_z,a_z,b_z,d_z_half,K_z_half,alpha_prime_z_half,a_z_half,b_z_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
+ double precision, dimension(NZ) :: d_z,K_z,alpha_z,a_z,b_z,d_z_half,K_z_half,alpha_z_half,a_z_half,b_z_half
! PML
double precision thickness_PML_x,thickness_PML_y,thickness_PML_z
@@ -419,8 +419,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -428,8 +428,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -437,8 +437,8 @@
d_z_half(:) = ZERO
K_z(:) = 1.d0
K_z_half(:) = 1.d0
- alpha_prime_z(:) = ZERO
- alpha_prime_z_half(:) = ZERO
+ alpha_z(:) = ZERO
+ alpha_z_half(:) = ZERO
a_z(:) = ZERO
a_z_half(:) = ZERO
@@ -463,7 +463,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -473,7 +473,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -488,7 +488,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -498,22 +498,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -538,7 +538,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -548,7 +548,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -563,7 +563,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -573,18 +573,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -609,7 +609,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -619,7 +619,7 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -634,7 +634,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -644,18 +644,18 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_prime_z(k)) * DELTAT)
- b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_prime_z_half(k)) * DELTAT)
+ b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_z(k)) * DELTAT)
+ b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_z_half(k)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_prime_z(k)))
+ if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_z(k)))
if(abs(d_z_half(k)) > 1.d-6) a_z_half(k) = d_z_half(k) * &
- (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_prime_z_half(k)))
+ (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_z_half(k)))
enddo
Modified: seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90 2011-05-28 13:54:47 UTC (rev 18483)
+++ seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90 2011-05-29 00:10:34 UTC (rev 18484)
@@ -212,8 +212,8 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 7.d0 ! from Gedney page 8.11
-! double precision, parameter :: ALPHA_MAX_PML = 0.d0 ! from festa and Vilotte
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+! double precision, parameter :: ALPHA_MAX_PML = 0.d0 ! from Festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -259,9 +259,9 @@
double precision :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz,div
! 1D arrays for the damping profiles
- double precision, dimension(1:NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(1:NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
- double precision, dimension(1:NZ) :: d_z,K_z,alpha_prime_z,a_z,b_z,d_z_half,K_z_half,alpha_prime_z_half,a_z_half,b_z_half
+ double precision, dimension(1:NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(1:NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
+ double precision, dimension(1:NZ) :: d_z,K_z,alpha_z,a_z,b_z,d_z_half,K_z_half,alpha_z_half,a_z_half,b_z_half
! PML
double precision thickness_PML_x,thickness_PML_y,thickness_PML_z
@@ -499,8 +499,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -508,8 +508,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -517,8 +517,8 @@
d_z_half(:) = ZERO
K_z(:) = 1.d0
K_z_half(:) = 1.d0
- alpha_prime_z(:) = ZERO
- alpha_prime_z_half(:) = ZERO
+ alpha_z(:) = ZERO
+ alpha_z_half(:) = ZERO
a_z(:) = ZERO
a_z_half(:) = ZERO
@@ -543,7 +543,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -553,7 +553,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -568,7 +568,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -578,22 +578,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -618,7 +618,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -628,7 +628,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -643,7 +643,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -653,18 +653,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -689,7 +689,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -699,7 +699,7 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -714,7 +714,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -724,18 +724,18 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_prime_z(k)) * DELTAT)
- b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_prime_z_half(k)) * DELTAT)
+ b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_z(k)) * DELTAT)
+ b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_z_half(k)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_prime_z(k)))
+ if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_z(k)))
if(abs(d_z_half(k)) > 1.d-6) a_z_half(k) = d_z_half(k) * &
- (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_prime_z_half(k)))
+ (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_z_half(k)))
enddo
More information about the CIG-COMMITS
mailing list