[cig-commits] [commit] devel: Zhinan fixed the real(kind=CUSTOM_REAL) vs double precision problems that remained in his CPML code and made the CPML code become unstable in some cases (cdc6cb7)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Apr 8 03:01:46 PDT 2014
Repository : ssh://geoshell/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/e4fa9d03bf2b0fc1837c42aa51eeb63f360575fe...fc67e6fd7ad890705b2b72b4b3c509accb22249e
>---------------------------------------------------------------
commit cdc6cb721d5ee9b2b6263de05cdb7e2e422fa04e
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Sat Mar 15 12:13:42 2014 +0100
Zhinan fixed the real(kind=CUSTOM_REAL) vs double precision problems that remained in his CPML code and made the CPML code become unstable in some cases
>---------------------------------------------------------------
cdc6cb721d5ee9b2b6263de05cdb7e2e422fa04e
src/specfem2D/compute_coupling_acoustic_el.f90 | 5 +--
src/specfem2D/compute_forces_acoustic.f90 | 4 +--
src/specfem2D/compute_forces_viscoelastic.F90 | 34 +++++++++---------
src/specfem2D/enforce_acoustic_free_surface.f90 | 6 ++--
src/specfem2D/invert_mass_matrix.F90 | 2 +-
src/specfem2D/pml_init.F90 | 47 +++++++------------------
src/specfem2D/specfem2D.F90 | 2 +-
7 files changed, 41 insertions(+), 59 deletions(-)
diff --git a/src/specfem2D/compute_coupling_acoustic_el.f90 b/src/specfem2D/compute_coupling_acoustic_el.f90
index 282243d..6e0a4e3 100644
--- a/src/specfem2D/compute_coupling_acoustic_el.f90
+++ b/src/specfem2D/compute_coupling_acoustic_el.f90
@@ -65,7 +65,7 @@
integer, dimension(num_fluid_solid_edges) :: fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge, &
fluid_solid_elastic_ispec,fluid_solid_elastic_iedge
integer :: nspec_PML
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
+ double precision, dimension(NGLLX,NGLLZ,nspec_PML) :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
logical:: PML_BOUNDARY_CONDITIONS
logical, dimension(nspec) :: is_PML
@@ -79,7 +79,8 @@
ispec_PML,CPML_region_local,singularity_type_xz
real(kind=CUSTOM_REAL) :: displ_x,displ_z,displ_n,&
xxi,zxi,xgamma,zgamma,jacobian1D,nx,nz,weight
- real(kind=CUSTOM_REAL) :: time,deltat,kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z,&
+ double precision :: time,deltat
+ double precision :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z, &
A8,A9,A10,bb_xz_1,bb_xz_2,coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2
! loop on all the coupling edges
diff --git a/src/specfem2D/compute_forces_acoustic.f90 b/src/specfem2D/compute_forces_acoustic.f90
index a2ab4f0..9452605 100644
--- a/src/specfem2D/compute_forces_acoustic.f90
+++ b/src/specfem2D/compute_forces_acoustic.f90
@@ -119,7 +119,7 @@
real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLZ,nspec_PML) :: rmemory_potential_acoustic
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML,2) :: &
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
+ double precision, dimension(NGLLX,NGLLZ,nspec_PML) :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
logical :: PML_BOUNDARY_CONDITIONS,STACEY_BOUNDARY_CONDITIONS
@@ -152,7 +152,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: potential_dot_dot_acoustic_PML
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: PML_dux_dxl,PML_dux_dzl,PML_dux_dxl_old,PML_dux_dzl_old
- real(kind=CUSTOM_REAL) :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z,time_n,time_nsub1,&
+ double precision :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z,time_n,time_nsub1,&
A5,A6,A7, bb_zx_1,bb_zx_2,coef0_zx_1,coef1_zx_1,coef2_zx_1,coef0_zx_2,coef1_zx_2,coef2_zx_2,&
A8,A9,A10,bb_xz_1,bb_xz_2,coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2,&
A0,A1,A2,A3,A4,bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2
diff --git a/src/specfem2D/compute_forces_viscoelastic.F90 b/src/specfem2D/compute_forces_viscoelastic.F90
index 8e26c15..acfd66b 100644
--- a/src/specfem2D/compute_forces_viscoelastic.F90
+++ b/src/specfem2D/compute_forces_viscoelastic.F90
@@ -135,7 +135,7 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) ::dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
!nsub1 denote discrete time step n-1
dux_dxl_nsub1,duz_dzl_nsub1,duz_dxl_nsub1,dux_dzl_nsub1
- real(kind=CUSTOM_REAL) :: coef0,coef1,coef2
+ double precision :: coef0,coef1,coef2
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
@@ -178,7 +178,8 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplusmul_relaxed_viscoel
! for attenuation
- real(kind=CUSTOM_REAL) :: phinu1,phinu2,tauinvnu1,tauinvnu2,theta_n_u,theta_nsub1_u
+ real(kind=CUSTOM_REAL) :: phinu1,phinu2,theta_n_u,theta_nsub1_u
+ double precision :: tauinvnu1,tauinvnu2
! for anisotropy
double precision :: c11,c15,c13,c33,c35,c55,c12,c23,c25
@@ -214,7 +215,7 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML,2) :: &
rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
+ double precision, dimension(NGLLX,NGLLZ,nspec_PML) :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ) :: accel_elastic_PML
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) ::PML_dux_dxl,PML_dux_dzl,PML_duz_dxl,PML_duz_dzl,&
@@ -223,7 +224,8 @@ subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
real(kind=CUSTOM_REAL) :: dux_dxi_old,dux_dgamma_old,duz_dxi_old,duz_dgamma_old
logical :: backward_simulation
- real(kind=CUSTOM_REAL) :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z,time_n,time_nsub1, &
+ double precision :: time_n,time_nsub1
+ double precision :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z, &
A5,A6,A7, bb_zx_1,bb_zx_2,coef0_zx_1,coef1_zx_1,coef2_zx_1,coef0_zx_2,coef1_zx_2,coef2_zx_2,&
A8,A9,A10,bb_xz_1,bb_xz_2,coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2,&
A0,A1,A2,A3,A4,bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2
@@ -1639,7 +1641,7 @@ end subroutine compute_forces_viscoelastic
implicit none
include "constants.h"
- real(kind=CUSTOM_REAL) :: bb,coef0,coef1,coef2
+ double precision :: bb,coef0,coef1,coef2
double precision :: deltat
coef0 = exp(- bb * deltat)
@@ -1662,17 +1664,17 @@ end subroutine compute_forces_viscoelastic
implicit none
include "constants.h"
- real(kind=CUSTOM_REAL), intent(in) :: time
+ double precision, intent(in) :: time
double precision :: deltat
- real(kind=CUSTOM_REAL), intent(in) :: kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z
+ double precision, intent(in) :: kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z
integer, intent(in) :: CPML_region_local,index_ik
- real(kind=CUSTOM_REAL), intent(out) :: A_0,A_1,A_2
- real(kind=CUSTOM_REAL), intent(out) :: coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2
+ double precision, intent(out) :: A_0,A_1,A_2
+ double precision, intent(out) :: coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2
integer, intent(out) :: singularity_type_2
!local variable
- real(kind=CUSTOM_REAL) :: bar_A_0,bar_A_1,bar_A_2,alpha_0,bb_1,bb_2
+ double precision :: bar_A_0,bar_A_1,bar_A_2,alpha_0,bb_1,bb_2
integer :: CPML_X_ONLY_TEMP,CPML_Z_ONLY_TEMP,CPML_XZ_ONLY_TEMP
logical,parameter :: FIRST_ORDER_CONVOLUTION = .false.
@@ -1765,18 +1767,18 @@ end subroutine compute_forces_viscoelastic
implicit none
include "constants.h"
- real(kind=CUSTOM_REAL), intent(in) :: time
+ double precision :: time
double precision :: deltat
- real(kind=CUSTOM_REAL), intent(in) :: kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z
+ double precision, intent(in) :: kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z
integer, intent(in) :: CPML_region_local
- real(kind=CUSTOM_REAL), intent(out) :: A_0, A_1, A_2, A_3, A_4
- real(kind=CUSTOM_REAL), intent(out) :: bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2
+ double precision, intent(out) :: A_0, A_1, A_2, A_3, A_4
+ double precision, intent(out) :: bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2
integer, intent(out) :: singularity_type
!local variable
- real(kind=CUSTOM_REAL) :: bar_A_0, bar_A_1, bar_A_2, bar_A_3, bar_A_4
- real(kind=CUSTOM_REAL) :: alpha_0, beta_xyz_1, beta_xyz_2
+ double precision :: bar_A_0, bar_A_1, bar_A_2, bar_A_3, bar_A_4
+ double precision :: alpha_0, beta_xyz_1, beta_xyz_2
beta_xyz_1 = beta_x + beta_z
beta_xyz_2 = beta_x * beta_z
diff --git a/src/specfem2D/enforce_acoustic_free_surface.f90 b/src/specfem2D/enforce_acoustic_free_surface.f90
index c24950e..cb8ee7d 100644
--- a/src/specfem2D/enforce_acoustic_free_surface.f90
+++ b/src/specfem2D/enforce_acoustic_free_surface.f90
@@ -80,9 +80,9 @@
iglob = ibool(i,j,ispec)
! make sure that an acoustic free surface is not enforced on periodic edges
if(.not. this_ibool_is_a_periodic_edge(iglob)) then
- potential_acoustic(iglob) = 0._CUSTOM_REAL
- potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
- potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ potential_acoustic(iglob) = ZERO
+ potential_dot_acoustic(iglob) = ZERO
+ potential_dot_dot_acoustic(iglob) = ZERO
endif
enddo
enddo
diff --git a/src/specfem2D/invert_mass_matrix.F90 b/src/specfem2D/invert_mass_matrix.F90
index 213b717..a74a04e 100644
--- a/src/specfem2D/invert_mass_matrix.F90
+++ b/src/specfem2D/invert_mass_matrix.F90
@@ -117,7 +117,7 @@
integer :: nspec_PML,ispec_PML
integer, dimension(nspec) :: spec_to_PML
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: K_x_store,K_z_store,d_x_store,d_z_store
+ double precision, dimension(NGLLX,NGLLZ,nspec_PML) :: K_x_store,K_z_store,d_x_store,d_z_store
logical, dimension(nspec) :: is_PML
logical :: PML_BOUNDARY_CONDITIONS,this_element_has_PML
integer, dimension(nspec) :: region_CPML
diff --git a/src/specfem2D/pml_init.F90 b/src/specfem2D/pml_init.F90
index 69580ba..c443adc 100644
--- a/src/specfem2D/pml_init.F90
+++ b/src/specfem2D/pml_init.F90
@@ -454,7 +454,7 @@ end subroutine pml_init
logical, dimension(nspec) :: is_PML
integer, dimension(nspec) :: region_CPML,spec_to_PML
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: K_x_store,K_z_store,d_x_store,d_z_store,&
+ double precision, dimension(NGLLX,NGLLZ,nspec_PML) :: K_x_store,K_z_store,d_x_store,d_z_store,&
alpha_x_store,alpha_z_store
! PML fixed parameters to compute parameter in PML
@@ -631,24 +631,12 @@ end subroutine pml_init
d0_z_bottom = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_z_bottom)
d0_z_top = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_z_top)
-! if (myrank == 0) then
-! write(IOUT,*)
-! write(IOUT,*) 'PML properties -------',myrank,'myrank'
-! write(IOUT,*) ' Vpmax=', vpmax
-! write(IOUT,*) ' log(Rcoef)=',log(Rcoef)
-! write(IOUT,*) ' thickness_PML_z_bottom =',thickness_PML_z_bottom
-! write(IOUT,*) ' thickness_PML_x_right =',thickness_PML_x_right
-! write(IOUT,*) ' thickness_PML_z_top =',thickness_PML_z_top
-! write(IOUT,*) ' thickness_PML_x_left =',thickness_PML_x_left
-! write(IOUT,*) ' d0_bottom =', d0_z_bottom
-! write(IOUT,*) ' d0_right =', d0_x_right
-! write(IOUT,*) ' d0_top =', d0_z_top
-! write(IOUT,*) ' d0_left =', d0_x_left
-! endif
-
- d_x_store = 0._CUSTOM_REAL; d_z_store = 0._CUSTOM_REAL
- K_x_store = 1._CUSTOM_REAL; K_z_store = 1._CUSTOM_REAL
- alpha_x_store = 0._CUSTOM_REAL; alpha_z_store = 0._CUSTOM_REAL
+ d_x_store = 0.d0
+ d_z_store = 0.d0
+ K_x_store = 1.d0
+ K_z_store = 1.d0
+ alpha_x_store = 0.d0
+ alpha_z_store = 0.d0
! define damping profile at the grid points
do ispec = 1,nspec
@@ -745,21 +733,12 @@ end subroutine pml_init
endif
endif
- if(CUSTOM_REAL == SIZE_REAL) then
- d_x_store(i,j,ispec_PML) = sngl(d_x)
- d_z_store(i,j,ispec_PML) = sngl(d_z)
- K_x_store(i,j,ispec_PML) = sngl(K_x)
- K_z_store(i,j,ispec_PML) = sngl(K_z)
- alpha_x_store(i,j,ispec_PML) = sngl(alpha_x)
- alpha_z_store(i,j,ispec_PML) = sngl(alpha_z)
- else
- d_x_store(i,j,ispec_PML) = d_x
- d_z_store(i,j,ispec_PML) = d_z
- K_x_store(i,j,ispec_PML) = K_x
- K_z_store(i,j,ispec_PML) = K_z
- alpha_x_store(i,j,ispec_PML) = alpha_x
- alpha_z_store(i,j,ispec_PML) = alpha_z
- endif
+ d_x_store(i,j,ispec_PML) = d_x
+ d_z_store(i,j,ispec_PML) = d_z
+ K_x_store(i,j,ispec_PML) = K_x
+ K_z_store(i,j,ispec_PML) = K_z
+ alpha_x_store(i,j,ispec_PML) = alpha_x
+ alpha_z_store(i,j,ispec_PML) = alpha_z
enddo; enddo
endif
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index ca41285..b90ccd6 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -998,7 +998,7 @@
!PML parameters
logical, dimension(:), allocatable :: is_PML
integer, dimension(:), allocatable :: region_CPML
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ double precision, dimension(:,:,:), allocatable :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
More information about the CIG-COMMITS
mailing list