Skip to content

Commit 258c3fc

Browse files
danielpeterdanielpeter
committed
moves PML boundary constants to setup/constants.h.in file
1 parent 5767c7b commit 258c3fc

13 files changed

+206
-127
lines changed

‎setup/constants.h.in‎

Lines changed: 22 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -230,12 +230,33 @@
230230
!!-----------------------------------------------------------
231231

232232
! power (exponent) used for the PML damping profile
233-
double precision, parameter :: NPOWER = 1.d0
233+
double precision, parameter :: CPML_NPOWER = 1.d0
234234

235235
! C-PML theoretical reflection coefficient
236236
! (INRIA research report section 6.1: http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf)
237237
double precision, parameter :: CPML_Rcoef = 0.001d0
238238

239+
! C-PML scaling factors K
240+
! the SPECFEM implementation uses
241+
! K_x = K_MIN + (K_MAX - ONE) * dist
242+
! Xie et al. (2014) paper formulates
243+
! K_x = K0 + K1 * x/L = K0 + K1 * dist
244+
! Thus, K0 == K_MIN and K1 == K_MAX - 1 or vice versa K_MAX == K1 + 1
245+
!
246+
! As mentioned in Xie et al. (2014, page 1729), as a default and for 2D cases they use
247+
! K0 == 1 and K1 == 0 -> that is equal to K_MIN == 1 and K_MAX == 1
248+
! and for distorted 3D meshes they use
249+
! K0 == 4 and K1 == 4 -> that is equal to K_MIN == 4 and K_MAX == 5
250+
double precision, parameter :: CPML_K_MIN = 1.d0 ! see also Gedney page 8.11
251+
double precision, parameter :: CPML_K_MAX = 1.d0
252+
253+
! C-PML convolution
254+
! by default, a second-order convolution scheme is used, turn this flag on for first-order scheme
255+
logical,parameter :: CPML_FIRST_ORDER_CONVOLUTION = .false.
256+
257+
! time scheme factor (see Xie et al. (2014), p. 1728, suggests 0 <= theta <= 1/4)
258+
double precision, parameter :: CPML_THETA = 1.d0/8.d0
259+
239260
! flags for the seven CPML regions
240261
integer, parameter :: CPML_X_ONLY = 1
241262
integer, parameter :: CPML_Y_ONLY = 2

‎src/generate_databases/generate_databases_par.F90‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -40,10 +40,7 @@ module generate_databases_par
4040
IMODEL_1D_PREM,IMODEL_1D_PREM_PB,IMODEL_1D_CASCADIA,IMODEL_1D_SOCAL, &
4141
IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC, &
4242
NX_TOPO_FILE,NY_TOPO_FILE, &
43-
ATTENUATION_COMP_MAXIMUM, &
44-
INJECTION_TECHNIQUE_IS_FK, &
45-
CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, &
46-
CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ,NPOWER,CPML_Rcoef
43+
ATTENUATION_COMP_MAXIMUM
4744

4845
use shared_parameters
4946

‎src/generate_databases/pml_set_local_dampingcoeff.f90‎

Lines changed: 128 additions & 79 deletions
Large diffs are not rendered by default.

‎src/gpu/kernels/UpdateDispVeloc_kernel.cu‎

Lines changed: 8 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,8 @@ __global__ void UpdateDispVeloc_PML_kernel(realw* displ,
7171
const int* d_CPML_to_spec,
7272
const int* d_ibool,
7373
const realw deltatsqover2,
74-
const realw deltatover2) {
74+
const realw deltatover2,
75+
const realw pml_theta) {
7576

7677
int ispec_cpml = blockIdx.x + blockIdx.y*gridDim.x;
7778
int ijk = threadIdx.x;
@@ -89,19 +90,17 @@ __global__ void UpdateDispVeloc_PML_kernel(realw* displ,
8990

9091
int iglob = d_ibool[ijk + NGLL3_PADDED*ispec] - 1;
9192

92-
const realw theta = 0.125f; // theta = 1.0 / 8.0;
93-
9493
// updates PML displacement
9594
PML_displ[INDEX5(NDIM,NGLLX,NGLLX,NGLLX,0,I,J,K,ispec_cpml)] = displ[3*iglob]
96-
+ deltatover2 * (1.0f - 2.0f * theta) * veloc[3*iglob]
97-
+ deltatsqover2 * (1.0f - theta) * accel[3*iglob];
95+
+ deltatover2 * (1.0f - 2.0f * pml_theta) * veloc[3*iglob]
96+
+ deltatsqover2 * (1.0f - pml_theta) * accel[3*iglob];
9897

9998
PML_displ[INDEX5(NDIM,NGLLX,NGLLX,NGLLX,1,I,J,K,ispec_cpml)] = displ[3*iglob+1]
100-
+ deltatover2 * (1.0f - 2.0f * theta) * veloc[3*iglob+1]
101-
+ deltatsqover2 * (1.0f - theta) * accel[3*iglob+1];
99+
+ deltatover2 * (1.0f - 2.0f * pml_theta) * veloc[3*iglob+1]
100+
+ deltatsqover2 * (1.0f - pml_theta) * accel[3*iglob+1];
102101

103102
PML_displ[INDEX5(NDIM,NGLLX,NGLLX,NGLLX,2,I,J,K,ispec_cpml)] = displ[3*iglob+2]
104-
+ deltatover2 * (1.0f - 2.0f * theta) * veloc[3*iglob+2]
105-
+ deltatsqover2 * (1.0f - theta) * accel[3*iglob+2];
103+
+ deltatover2 * (1.0f - 2.0f * pml_theta) * veloc[3*iglob+2]
104+
+ deltatsqover2 * (1.0f - pml_theta) * accel[3*iglob+2];
106105
}
107106
}

‎src/gpu/kernels/kernel_proto.cu.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1305,7 +1305,8 @@ __global__ void UpdateDispVeloc_PML_kernel(realw* displ,
13051305
const int* d_CPML_to_spec,
13061306
const int* d_ibool,
13071307
const realw deltatsqover2,
1308-
const realw deltatover2) ;
1308+
const realw deltatover2,
1309+
const realw pml_theta) ;
13091310

13101311

13111312
//

‎src/gpu/mesh_constants_gpu.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -737,6 +737,8 @@ typedef struct mesh_ {
737737
realw* d_pml_convolution_coef_strain;
738738
realw* d_pml_convolution_coef_abar;
739739

740+
realw pml_theta; // passed from constants.h (by default CPML_THETA = 1.0 / 8.0)
741+
740742
// ------------------------------------------------------------------ //
741743
// acoustic wavefield
742744
// ------------------------------------------------------------------ //

‎src/gpu/prepare_mesh_constants_cuda.cu‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1202,7 +1202,8 @@ void FC_FUNC_(prepare_fields_elastic_pml,
12021202
realw* pml_convolution_coef_abar,
12031203
realw* pml_convolution_coef_strain,
12041204
realw* h_wgll_cube,
1205-
realw* rhostore) {
1205+
realw* rhostore,
1206+
double* CPML_THETA) {
12061207

12071208
TRACE("prepare_fields_elastic_pml");
12081209
int size;
@@ -1212,8 +1213,12 @@ void FC_FUNC_(prepare_fields_elastic_pml,
12121213
// checks if anything to do
12131214
if (! mp->pml_conditions) return;
12141215

1216+
// PML elements
12151217
mp->NSPEC_CPML = *NSPEC_CPML;
12161218

1219+
// PML time scheme factor (passed as argument to be consistent with chosen setting in constants.h)
1220+
mp->pml_theta = (realw) *CPML_THETA;
1221+
12171222
gpuCreateCopy_todevice_int((void**)&mp->d_is_CPML,is_CPML,mp->NSPEC_AB);
12181223

12191224
if (mp->NSPEC_CPML > 0){

‎src/gpu/specfem3D_gpu_cuda_method_stubs.c‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -707,7 +707,8 @@ void FC_FUNC_(prepare_fields_elastic_pml,
707707
realw* pml_convolution_coef_abar,
708708
realw* pml_convolution_coef_strain,
709709
realw* h_wgll_cube,
710-
realw* rhostore) {}
710+
realw* rhostore,
711+
double* CPML_THETA) {}
711712

712713
void FC_FUNC_(prepare_sim2_or_3_const_device,
713714
PREPARE_SIM2_OR_3_CONST_DEVICE)(long* Mesh_pointer,int *nadj_rec_local, int* NTSTEP_BETWEEN_READ_ADJSRC,
@@ -1048,10 +1049,10 @@ void FC_FUNC_(unregister_host_array,
10481049

10491050
void FC_FUNC_(update_displacement_cuda,
10501051
UPDATE_DISPLACMENT_CUDA)(long* Mesh_pointer,
1051-
realw* deltat_F,
1052-
realw* deltatover2_F,
1053-
realw* deltatsqover2_F,
1054-
int* FORWARD_OR_ADJOINT) {}
1052+
realw* deltat_F,
1053+
realw* deltatover2_F,
1054+
realw* deltatsqover2_F,
1055+
int* FORWARD_OR_ADJOINT) {}
10551056

10561057
void FC_FUNC_(update_displacement_ac_cuda,
10571058
UPDATE_DISPLACEMENT_AC_CUDA)(long* Mesh_pointer,

‎src/gpu/update_displacement_cuda.cu‎

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -38,10 +38,10 @@
3838
extern EXTERN_LANG
3939
void FC_FUNC_(update_displacement_cuda,
4040
UPDATE_DISPLACMENT_CUDA)(long* Mesh_pointer,
41-
realw* deltat_F,
42-
realw* deltatover2_F,
43-
realw* deltatsqover2_F,
44-
int* FORWARD_OR_ADJOINT) {
41+
realw* deltat_F,
42+
realw* deltatover2_F,
43+
realw* deltatsqover2_F,
44+
int* FORWARD_OR_ADJOINT) {
4545

4646
TRACE("update_displacement_cuda");
4747

@@ -108,7 +108,8 @@ void FC_FUNC_(update_displacement_cuda,
108108
mp->d_PML_displ_old,
109109
mp->NSPEC_CPML,mp->d_CPML_to_spec,
110110
mp->d_ibool,
111-
deltatsqover2,deltatover2);
111+
deltatsqover2,deltatover2,
112+
mp->pml_theta);
112113
}
113114
#endif
114115
#ifdef USE_HIP
@@ -118,7 +119,8 @@ void FC_FUNC_(update_displacement_cuda,
118119
mp->d_PML_displ_old,
119120
mp->NSPEC_CPML,mp->d_CPML_to_spec,
120121
mp->d_ibool,
121-
deltatsqover2,deltatover2);
122+
deltatsqover2,deltatover2,
123+
mp->pml_theta);
122124
}
123125
#endif
124126
}
@@ -149,7 +151,8 @@ void FC_FUNC_(update_displacement_cuda,
149151
mp->d_PML_displ_new,
150152
mp->NSPEC_CPML,mp->d_CPML_to_spec,
151153
mp->d_ibool,
152-
deltatsqover2,deltatover2);
154+
deltatsqover2,deltatover2,
155+
mp->pml_theta);
153156
}
154157
#endif
155158
#ifdef USE_HIP
@@ -159,7 +162,8 @@ void FC_FUNC_(update_displacement_cuda,
159162
mp->d_PML_displ_new,
160163
mp->NSPEC_CPML,mp->d_CPML_to_spec,
161164
mp->d_ibool,
162-
deltatsqover2,deltatover2);
165+
deltatsqover2,deltatover2,
166+
mp->pml_theta);
163167
}
164168
#endif
165169

‎src/specfem3D/pml_par.f90‎

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ module pml_par
3030

3131
! main parameter module for C-PML simulations
3232

33-
use constants, only: CUSTOM_REAL
33+
use constants, only: CUSTOM_REAL,CPML_THETA
3434

3535
implicit none
3636

@@ -69,13 +69,15 @@ module pml_par
6969
! minimum distance between parameters of CPML to avoid the singularities
7070
real(kind=CUSTOM_REAL) :: min_distance_between_CPML_parameter
7171

72+
! PML time scheme factors
73+
real(kind=CUSTOM_REAL), parameter :: ONE_MINUS_THETA = real(1.d0 - CPML_THETA,kind=CUSTOM_REAL)
74+
real(kind=CUSTOM_REAL), parameter :: ONE_MINUS_TWO_THETA = real(1.d0 - 2.d0 * CPML_THETA,kind=CUSTOM_REAL)
75+
7276
! store the field of displ + (1-2 * \theta)/2*deltat * veloc + (1-\theta)/2*deltat**2 * accel
7377
! for second order convolution scheme
7478
! where displ, veloc, accel are defined at n-1 time step
7579
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: PML_displ_old
7680

77-
real(kind=CUSTOM_REAL), parameter :: THETA = real(1.d0/8.d0, kind=CUSTOM_REAL)
78-
7981
! store the field of displ + (1-2 * \theta)/2*deltat * veloc for second order convolution scheme
8082
! where displ is defined at n time step, while veloc is predicted veloc at "n" time step
8183
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: PML_displ_new
@@ -134,7 +136,4 @@ module pml_par
134136
integer :: b_reclen_PML_field,b_reclen_PML_potential
135137
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_PML_field,b_PML_potential
136138

137-
! convolution coefficients
138-
logical,parameter :: FIRST_ORDER_CONVOLUTION = .false.
139-
140139
end module pml_par

0 commit comments

Comments
 (0)