Skip to content

Commit ac7e451

Browse files
danielpeterdanielpeter
committed
adds acoustic-elastic coupling for PML elements on GPU
1 parent 17b65b9 commit ac7e451

10 files changed

Lines changed: 522 additions & 211 deletions

‎src/gpu/compute_coupling_cuda.cu‎

Lines changed: 26 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -105,7 +105,18 @@ void FC_FUNC_(compute_coupling_ac_el_cuda,
105105
mp->d_coupling_ac_el_jacobian2Dw,
106106
mp->d_ibool,
107107
mp->simulation_type,
108-
backward_simulation);
108+
backward_simulation,
109+
mp->pml_boundary_conditions,
110+
mp->d_spec_to_pml,
111+
mp->nspec_pml_x,
112+
mp->nspec_pml_z,
113+
mp->deltat,
114+
mp->d_displ_elastic_old,
115+
mp->d_rmemory_fsb_displ_elastic,
116+
mp->alphax_store,
117+
mp->alphaz_store,
118+
mp->betax_store,
119+
mp->betaz_store);
109120

110121
//double end_time = get_time();
111122
//printf("Elapsed time: %e\n",end_time-start_time);
@@ -183,7 +194,20 @@ void FC_FUNC_(compute_coupling_el_ac_cuda,
183194
mp->d_coupling_ac_el_jacobian2Dw,
184195
mp->d_ibool,
185196
mp->simulation_type,
186-
backward_simulation);
197+
backward_simulation,
198+
mp->pml_boundary_conditions,
199+
mp->d_spec_to_pml,
200+
mp->nspec_pml_x,
201+
mp->nspec_pml_z,
202+
mp->deltat,
203+
mp->d_potential_acoustic,
204+
mp->d_potential_dot_acoustic,
205+
mp->d_potential_old,
206+
mp->d_rmemory_sfb_potential_ddot_acoustic,
207+
mp->alphax_store,
208+
mp->alphaz_store,
209+
mp->betax_store,
210+
mp->betaz_store);
187211

188212
//double end_time = get_time();
189213
//printf("Elapsed time: %e\n",end_time-start_time);

‎src/gpu/kernels/compute_coupling_acoustic_el_kernel.cu‎

Lines changed: 121 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,18 @@ __global__ void compute_coupling_acoustic_el_kernel(realw* displ,
4343
realw* coupling_ac_el_jacobian1Dw,
4444
int* d_ibool,
4545
int simulation_type,
46-
int backward_simulation) {
46+
int backward_simulation,
47+
const int PML,
48+
const int* d_spec_to_pml,
49+
const int NSPEC_PML_X,
50+
const int NSPEC_PML_Z,
51+
const realw deltat,
52+
realw_const_p d_displ_elastic_old,
53+
realw* d_rmemory_fsb_displ_elastic,
54+
realw_const_p alphax_store,
55+
realw_const_p alphaz_store,
56+
realw_const_p betax_store,
57+
realw_const_p betaz_store) {
4758

4859
int igll = threadIdx.x;
4960
int iface = blockIdx.x + gridDim.x*blockIdx.y;
@@ -54,11 +65,6 @@ __global__ void compute_coupling_acoustic_el_kernel(realw* displ,
5465
realw jacobianw;
5566

5667
if (iface < num_coupling_ac_el_faces){
57-
58-
// don't compute points outside NGLLSQUARE==NGLL2==25
59-
// way 2: no further check needed since blocksize = 25
60-
// if (igll<NGLL2) {
61-
6268
// "-1" from index values to convert from Fortran-> C indexing
6369
ispec = coupling_ac_el_ispec[iface] - 1;
6470

@@ -79,13 +85,120 @@ __global__ void compute_coupling_acoustic_el_kernel(realw* displ,
7985
displ_z = - displ_z;
8086
}
8187

88+
// PML
89+
if (PML) {
90+
// PML element index
91+
int ispec_pml = d_spec_to_pml[ispec] - 1;
92+
// checks if element is inside the PML
93+
if (ispec_pml >= 0) {
94+
realw alpha1,beta1;
95+
realw coef0_1,coef1_1,coef2_1;
96+
realw A9;
97+
98+
// to match offset for local (i,j) index to thread index (tx == i + j * NGLLX from I = (tx-J*NGLLX))
99+
int tx = i + j * NGLLX;
100+
101+
// local PML array index
102+
int offset_pml = ispec_pml * NGLL2 + tx; // ispec_pml elements in range [0,NSPEC_PML-1]
103+
104+
// coefficients
105+
if (ispec_pml < NSPEC_PML_X){
106+
// in CPML_X_ONLY
107+
// alpha1 == alpha_x
108+
// alpha_z == 0
109+
//
110+
// beta1 == beta_x == alpha_x + d_x / K_x
111+
// beta_z == 0
112+
//
113+
alpha1 = alphax_store[offset_pml];
114+
beta1 = betax_store[offset_pml];
115+
} else if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
116+
// in CPML_Z_ONLY region
117+
// alpha1 == alpha_z
118+
// alpha_x == 0
119+
//
120+
// beta1 == beta_z == alpha_z + d_z / K_z
121+
// beta_x == 0
122+
alpha1 = alphaz_store[offset_pml];
123+
beta1 = betaz_store[offset_pml];
124+
} else {
125+
// not used, there should be no coupling interface in the CPML_XZ regions at the corners
126+
alpha1 = 0.f;
127+
beta1 = 0.f;
128+
}
129+
130+
// for all PML regions
131+
realw c1 = __expf(-0.5f * deltat * alpha1);
132+
133+
coef0_1 = c1 * c1;
134+
if (abs(alpha1) > 0.00001f){
135+
// coef1_zx_1 == (1 - c1)/alpha1
136+
// coef2_zx_1 == coef1 * c1
137+
coef1_1 = (1.f - c1) / alpha1;
138+
coef2_1 = coef1_1 * c1;
139+
} else {
140+
// coef1_zx_1 == 1/2 dt
141+
// coef2_zx_1 == coef1_zx_1
142+
coef1_1 = 0.5f * deltat;
143+
coef2_1 = coef1_1;
144+
}
145+
146+
// memory variables update
147+
// see compute_coupling_acoustic_el.f90 (line ~142):
148+
// ! Newmark
149+
// rmemory_fsb_displ_elastic(1,1,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,1,i,j,inum) + &
150+
// coef1_xz_1 * displ_elastic(1,iglob) + coef2_xz_1 * displ_elastic_old(1,iglob)
151+
// rmemory_fsb_displ_elastic(1,2,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,2,i,j,inum) + &
152+
// coef1_xz_1 * displ_elastic(2,iglob) + coef2_xz_1 * displ_elastic_old(2,iglob)
153+
//
154+
// x-comp
155+
realw r_x = coef0_1 * d_rmemory_fsb_displ_elastic[INDEX3(NDIM,NGLLX,0,igll,iface)] + coef1_1 * displ_x + coef2_1 * d_displ_elastic_old[offset_pml*2];
156+
// z-comp
157+
realw r_z = coef0_1 * d_rmemory_fsb_displ_elastic[INDEX3(NDIM,NGLLX,1,igll,iface)] + coef1_1 * displ_z + coef2_1 * d_displ_elastic_old[offset_pml*2+1];
158+
159+
d_rmemory_fsb_displ_elastic[INDEX3(NDIM,NGLLX,0,igll,iface)] = r_x; // (1,igll,iface)
160+
d_rmemory_fsb_displ_elastic[INDEX3(NDIM,NGLLX,1,igll,iface)] = r_z; // (2,igll,iface)
161+
162+
// displacement update
163+
// see compute_coupling_acoustic_el.f90 (line ~166):
164+
// displ_x = A8 * displ_elastic(1,iglob) + A9 * rmemory_fsb_displ_elastic(1,1,i,j,inum)
165+
// displ_z = A8 * displ_elastic(2,iglob) + A9 * rmemory_fsb_displ_elastic(1,2,i,j,inum)
166+
//
167+
// with coefficients A8, A9 from routine lik_parameter_computation(..),
168+
// note that we require K_x == K_z == 1:
169+
// for CPML_X_ONLY_TEMP
170+
// A8 == A_0 == kappa_x == 1
171+
// A9 == - A_0 * (alpha_x - beta_x) == beta_x - alpha_x
172+
// for CPML_Z_ONLY_TEMP
173+
// A8 == A_0 == 1 / kappa_z == 1
174+
// A9 == - A_0 * (beta_z - alpha_z) == alpha_z - beta_z
175+
if (ispec_pml < NSPEC_PML_X){
176+
// in CPML_X_ONLY region
177+
//A8 = 1.0f;
178+
A9 = beta1 - alpha1;
179+
} else if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)) {
180+
// in CPML_Z_ONLY region
181+
//A8 = 1.0f;
182+
A9 = alpha1 - beta1;
183+
} else {
184+
// in CPML_XZ region
185+
// should not occur
186+
//A8 = 1.0f; // keeps displ_x and displ_z as is
187+
A9 = 0.f;
188+
}
189+
// overwrites displ_x and displ_z
190+
displ_x += A9 * r_x;
191+
displ_z += A9 * r_z;
192+
}
193+
} // PML
194+
82195
// gets associated normal on GLL point
83196
nx = coupling_ac_el_normal[INDEX3(NDIM,NGLLX,0,igll,iface)]; // (1,igll,iface)
84197
nz = coupling_ac_el_normal[INDEX3(NDIM,NGLLX,1,igll,iface)]; // (2,igll,iface)
85198

86199
// calculates displacement component along normal
87200
// (normal points outwards of acoustic element)
88-
displ_n = displ_x*nx + displ_z*nz;
201+
displ_n = displ_x * nx + displ_z * nz;
89202

90203
// gets associated, weighted jacobian
91204
jacobianw = coupling_ac_el_jacobian1Dw[INDEX2(NGLLX,igll,iface)];
@@ -99,9 +212,7 @@ __global__ void compute_coupling_acoustic_el_kernel(realw* displ,
99212
// (see e.g. Chaljub & Vilotte, Nissen-Meyer thesis...)
100213
// it also means you have to calculate and update this here first before
101214
// calculating the coupling on the elastic side for the acceleration...
102-
atomicAdd(&potential_dot_dot_acoustic[iglob],+ jacobianw*displ_n);
103-
104-
// }
215+
atomicAdd(&potential_dot_dot_acoustic[iglob],jacobianw * displ_n);
105216
}
106217
}
107218

‎src/gpu/kernels/compute_coupling_elastic_ac_kernel.cu‎

Lines changed: 137 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,20 @@ __global__ void compute_coupling_elastic_ac_kernel(realw* potential_dot_dot_acou
4343
realw* coupling_ac_el_jacobian1Dw,
4444
int* d_ibool,
4545
int simulation_type,
46-
int backward_simulation) {
46+
int backward_simulation,
47+
const int PML,
48+
const int* d_spec_to_pml,
49+
const int NSPEC_PML_X,
50+
const int NSPEC_PML_Z,
51+
const realw deltat,
52+
realw_const_p potential_acoustic,
53+
realw_const_p potential_dot_acoustic,
54+
realw_const_p d_potential_old,
55+
realw* d_rmemory_sfb_potential_ddot_acoustic,
56+
realw_const_p alphax_store,
57+
realw_const_p alphaz_store,
58+
realw_const_p betax_store,
59+
realw_const_p betaz_store) {
4760

4861
int igll = threadIdx.x;
4962
int iface = blockIdx.x + gridDim.x*blockIdx.y;
@@ -53,13 +66,7 @@ __global__ void compute_coupling_elastic_ac_kernel(realw* potential_dot_dot_acou
5366
realw nx,nz;
5467
realw jacobianw;
5568

56-
5769
if (iface < num_coupling_ac_el_faces){
58-
59-
// don't compute points outside NGLLSQUARE==NGLL2==25
60-
// way 2: no further check needed since blocksize = 25
61-
// if (igll<NGLL2) {
62-
6370
// "-1" from index values to convert from Fortran-> C indexing
6471
ispec = coupling_ac_el_ispec[iface] - 1;
6572

@@ -86,6 +93,127 @@ __global__ void compute_coupling_elastic_ac_kernel(realw* potential_dot_dot_acou
8693
pressure = - pressure;
8794
}
8895

96+
// PML
97+
if (PML) {
98+
// PML element index
99+
int ispec_pml = d_spec_to_pml[ispec] - 1;
100+
// checks if element is inside the PML
101+
if (ispec_pml >= 0) {
102+
realw alpha1,beta1;
103+
realw coef0_1,coef1_1,coef2_1;
104+
realw A1,A2,A3;
105+
106+
// to match offset for local (i,j) index to thread index (tx == i + j * NGLLX from I = (tx-J*NGLLX))
107+
int tx = i + j * NGLLX;
108+
109+
// local PML array index (tx == igll)
110+
int offset_pml = ispec_pml * NGLL2 + tx; // ispec_pml elements in range [0,NSPEC_PML-1]
111+
112+
// coefficients
113+
if (ispec_pml < NSPEC_PML_X){
114+
// in CPML_X_ONLY
115+
// alpha1 == alpha_x
116+
// alpha_z == 0
117+
//
118+
// beta1 == beta_x == alpha_x + d_x / K_x
119+
// beta_z == 0
120+
//
121+
alpha1 = alphax_store[offset_pml];
122+
beta1 = betax_store[offset_pml];
123+
} else if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
124+
// in CPML_Z_ONLY region
125+
// alpha1 == alpha_z
126+
// alpha_x == 0
127+
//
128+
// beta1 == beta_z == alpha_z + d_z / K_z
129+
// beta_x == 0
130+
alpha1 = alphaz_store[offset_pml];
131+
beta1 = betaz_store[offset_pml];
132+
} else {
133+
// should not occur - there should be no coupling edge interfaces in the CPML_XZ regions at the corners
134+
alpha1 = 0.f;
135+
beta1 = 0.f;
136+
}
137+
138+
// for all PML regions
139+
realw c1 = __expf(-0.5f * deltat * alpha1);
140+
141+
coef0_1 = c1 * c1;
142+
if (abs(alpha1) > 0.00001f){
143+
// coef1_zx_1 == (1 - c1)/alpha1
144+
// coef2_zx_1 == coef1 * c1
145+
coef1_1 = (1.f - c1) / alpha1;
146+
coef2_1 = coef1_1 * c1;
147+
} else {
148+
// coef1_zx_1 == 1/2 dt
149+
// coef2_zx_1 == coef1_zx_1
150+
coef1_1 = 0.5f * deltat;
151+
coef2_1 = coef1_1;
152+
}
153+
154+
// memory variables update
155+
// see compute_coupling_viscoelastic_ac.f90 (line ~127):
156+
// ! Newmark
157+
// rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) = &
158+
// coef0_1 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + &
159+
// coef1_1 * potential_acoustic(iglob) + coef2_1 * potential_acoustic_old(iglob)
160+
//
161+
realw rm = coef0_1 * d_rmemory_sfb_potential_ddot_acoustic[INDEX2(NGLLX,igll,iface)]
162+
+ coef1_1 * potential_acoustic[iglob] + coef2_1 * d_potential_old[offset_pml];
163+
164+
d_rmemory_sfb_potential_ddot_acoustic[INDEX2(NGLLX,igll,iface)] = rm;
165+
166+
// potential update
167+
// see compute_coupling_viscoelastic_ac.f90 (line ~144):
168+
// pressure = - (A0 * potential_dot_dot_acoustic(iglob) + A1 * potential_dot_acoustic(iglob) + &
169+
// A2 * potential_acoustic(iglob) + A3 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum))
170+
//
171+
// with coefficients A0,A1,A2,A3 from routine l_parameter_computation(..),
172+
// note that we require K_x == K_z == 1:
173+
// for CPML_X_ONLY
174+
// A_0 == kappa_x == 1
175+
// A_1 == A_0 * (beta_x - alpha_x)
176+
// A_2 == - A_0 * alpha_x * (beta_x - alpha_x)
177+
// A_3 == A_0 * alpha_x**2 * (beta_x - alpha_x)
178+
// A_4 == 0
179+
// for CPML_Z_ONLY_TEMP
180+
// A_0 == kappa_z == 1
181+
// A_1 == A_0 * (beta_z - alpha_z)
182+
// A_2 == - A_0 * alpha_z * (beta_z - alpha_z)
183+
// A_3 == 0
184+
// A_4 == A_0 * alpha_z**2 * (beta_z - alpha_z)
185+
if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
186+
// in CPML_X_ONLY or in CPML_Z_ONLY region
187+
//A0 = 1.0f;
188+
A1 = beta1 - alpha1;
189+
A2 = - alpha1 * A1; // - alpha1 * (beta1 - alpha1)
190+
A3 = - alpha1 * A2; // alpha1 * alpha1 * (beta1 - alpha1)
191+
//A4 = 0.f;
192+
193+
/*
194+
// note: there might be some difference with the CPU update of routine compute_coupling_viscoelastic_ac(), i.e,
195+
// for CPML_X_ONLY the coefficient A4 == 0 and for CPML_Z_ONLY the coefficient A3 == 0
196+
// thus, to match CPU implementation
197+
if (ispec_pml < NSPEC_PML_X){
198+
// in CPML_X_ONLY region
199+
A3 = - alpha1 * A2;
200+
} else {
201+
// in CPML_Z_ONLY region
202+
A3 = 0.f;
203+
}
204+
*/
205+
} else {
206+
// in CPML_XZ region - should not occur
207+
//A0 = 1.0f; // keeps pressure as is
208+
A1 = 0.f;
209+
A2 = 0.f;
210+
A3 = 0.f;
211+
}
212+
// overwrites pressure value
213+
pressure = - (potential_dot_dot_acoustic[iglob] + A1 * potential_dot_acoustic[iglob] + A2 * potential_acoustic[iglob] + A3 * rm);
214+
}
215+
} // PML
216+
89217
// continuity of displacement and pressure on global point
90218
//
91219
// note: Newmark time scheme together with definition of scalar potential:
@@ -95,10 +223,8 @@ __global__ void compute_coupling_elastic_ac_kernel(realw* potential_dot_dot_acou
95223
// (see e.g. Chaljub & Vilotte, Nissen-Meyer thesis...)
96224
// it means you have to calculate and update the acoustic pressure first before
97225
// calculating this term...
98-
atomicAdd(&accel[iglob*2],+ jacobianw*nx*pressure);
99-
atomicAdd(&accel[iglob*2+1],+ jacobianw*nz*pressure);
100-
101-
// }
226+
atomicAdd(&accel[iglob*2],jacobianw * nx * pressure);
227+
atomicAdd(&accel[iglob*2+1],jacobianw * nz * pressure);
102228
}
103229
}
104230

0 commit comments

Comments
 (0)