Skip to content

Commit ec2d201

Browse files
danielpeterdanielpeter
committed
adds PML on GPU for elastic elements
1 parent d9c2084 commit ec2d201

10 files changed

Lines changed: 562 additions & 102 deletions

‎src/gpu/compute_forces_viscoelastic_cuda.cu‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -304,6 +304,23 @@ TRACE("Kernel_2");
304304
mp->nspec_pml_x,
305305
mp->nspec_pml_z,
306306
mp->deltat,
307+
mp->PML_dux_dxl_old,
308+
mp->PML_dux_dzl_old,
309+
mp->PML_duz_dxl_old,
310+
mp->PML_duz_dzl_old,
311+
mp->d_displ_elastic_old,
312+
mp->d_rmemory_dux_dx,
313+
mp->d_rmemory_dux_dx2,
314+
mp->d_rmemory_duz_dx,
315+
mp->d_rmemory_duz_dx2,
316+
mp->d_rmemory_dux_dz,
317+
mp->d_rmemory_dux_dz2,
318+
mp->d_rmemory_duz_dz,
319+
mp->d_rmemory_duz_dz2,
320+
mp->d_rmemory_displ_elastic,
321+
mp->d_rmemory_displ_elastic2,
322+
mp->d_veloc,
323+
mp->d_rhostore,
307324
mp->alphax_store,
308325
mp->alphaz_store,
309326
mp->betax_store,

‎src/gpu/kernels/Kernel_2_acoustic_impl.cu‎

Lines changed: 12 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -463,7 +463,6 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
463463
rho_invl_times_jacobianl = 1.f /(rhol * (xixl*gammazl-gammaxl*xizl));
464464

465465
// loads hprime into shared memory
466-
467466
#ifdef USE_TEXTURES_CONSTANTS
468467
sh_hprime_xx[tx] = tex1Dfetch(d_hprime_xx_tex,tx);
469468
#else
@@ -767,13 +766,14 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
767766
// A5 == A0 == 1
768767
// A6 == A2 == - (beta_z - alpha_z) == (alpha_z - beta_z)
769768
// A7 == A1 == 0
770-
dpotentialdxl += (alpha1-beta1) * r1;
769+
realw bar_A = (alpha1-beta1);
770+
dpotentialdxl += bar_A * r1;
771771
// dux_dz: uses coefficients for X_ONLY_TEMP case
772772
// with arguments alpha_x -> alpha_x == alpha1, beta_x -> beta_x == beta1 (index_ik 13)
773773
// A8 == A0 == 1
774774
// A9 == A1 == - (alpha_x - beta_x)
775775
// A10 == A2 == 0
776-
dpotentialdzl -= (alpha1-beta1) * r2;
776+
dpotentialdzl -= bar_A * r2;
777777
} else if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)) {
778778
// in CPML_Z_ONLY region
779779
// (alpha1 == alpha_z and beta1 == beta_z)
@@ -783,13 +783,14 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
783783
// A5 == A0 == 1
784784
// A6 == A1 == - (alpha_x - beta_x)
785785
// A7 == A2 == 0
786-
dpotentialdxl -= (alpha1-beta1) * r1;
786+
realw bar_A = (alpha1-beta1);
787+
dpotentialdxl -= bar_A * r1;
787788
// dux_dz: uses coefficients for Z_ONLY_TEMP case
788789
// with arguments alpha_z -> alpha_z == alpha1, beta_z -> beta_z == beta1 (index_ik 13)
789790
// A8 == A0 == 1
790791
// A9 == A2 == - (beta_z - alpha_z) == (alpha_z - beta_z)
791792
// A10 == A1 == 0
792-
dpotentialdzl += (alpha1-beta1) * r2;
793+
dpotentialdzl += bar_A * r2;
793794
} else {
794795
// in CPML_XZ region
795796
// (alpha1 == alpha_z, beta1 == beta_z, and alphax == alpha_x, betax == beta_x)
@@ -800,24 +801,22 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
800801
// A5 == A0 == 1
801802
// A6 == A1 == 1/2 * (gamma_x - alpha_x)
802803
// gamma_x == (alpha_x * beta_z + alpha_x**2 + 2 * beta_x * alpha_z - 2 * alpha_x * (beta_x + alpha_z)) / (beta_z - alpha_x)
803-
dpotentialdxl += 0.5f * ((alpha1 * betax + alpha1*alpha1 + 2.f * beta1 * alphax - 2.f * alpha1 * (beta1 + alphax)) / (betax - alpha1)
804-
- alpha1) * r1;
804+
realw bar_A1 = 0.5f * ((alpha1 * betax + alpha1*alpha1 + 2.f * beta1 * alphax - 2.f * alpha1 * (beta1 + alphax)) / (betax - alpha1) - alpha1);
805805
// A7 == A2 == 1/2 * (gamma_z - beta_z)
806806
// gamma_z == (alpha_x * beta_z + beta_z**2 + 2 * beta_x * alpha_z - 2 * beta_z * (beta_x + alpha_z)) / (alpha_x - beta_z)
807-
dpotentialdxl += 0.5f * ((alpha1 * betax + betax*betax + 2.f * beta1 * alphax - 2.f * betax * ( beta1 + alphax)) / (alpha1 - betax)
808-
- betax) * r3;
807+
realw bar_A2 = 0.5f * ((alpha1 * betax + betax*betax + 2.f * beta1 * alphax - 2.f * betax * ( beta1 + alphax)) / (alpha1 - betax) - betax);
808+
dpotentialdxl += bar_A1 * r1 + bar_A2 * r3;
809809
// dux_dz: uses coefficients for XZ_TEMP case
810810
// with arguments alpha_x -> alpha_x == alphax, beta_x -> beta_x == betax (index_ik 13)
811811
// alpha_z -> alpha_z == alpha1, beta_z -> beta_z == beta1
812812
// A8 == A0 == 1
813813
// A9 == A2 == 1/2 * (gamma_z - beta_z)
814814
// gamma_z == (alpha_x * beta_z + beta_z**2 + 2 * beta_x * alpha_z - 2 * beta_z * (beta_x + alpha_z)) / (alpha_x - beta_z)
815-
dpotentialdzl += 0.5f * ((alphax * beta1 + beta1*beta1 + 2.f * betax * alpha1 - 2.f * beta1 * (betax + alpha1)) / (alphax - beta1)
816-
- beta1) * r2;
815+
realw bar_A3 = 0.5f * ((alphax * beta1 + beta1*beta1 + 2.f * betax * alpha1 - 2.f * beta1 * (betax + alpha1)) / (alphax - beta1) - beta1);
817816
// A10 == A1 == 1/2 * (gamma_x - alpha_x)
818817
// gamma_x == (alpha_x * beta_z + alpha_x**2 + 2 * beta_x * alpha_z - 2 * alpha_x * (beta_x + alpha_z)) / (beta_z - alpha_x)
819-
dpotentialdzl += 0.5f * ((alphax * beta1 + alphax*alphax + 2.f * betax * alpha1 - 2.f * alphax * (betax + alpha1)) / (beta1 - alphax)
820-
- alphax) * r4;
818+
realw bar_A4 = 0.5f * ((alphax * beta1 + alphax*alphax + 2.f * betax * alpha1 - 2.f * alphax * (betax + alpha1)) / (beta1 - alphax) - alphax);
819+
dpotentialdzl += bar_A3 * r2 + bar_A4 * r4;
821820
}
822821
//__syncthreads(); // not needed... still everything thread local
823822

0 commit comments

Comments
 (0)