Skip to content

Commit 17b65b9

Browse files
danielpeterdanielpeter
committed
updates usage of coefficient arrays for PML on GPU
1 parent ec2d201 commit 17b65b9

11 files changed

Lines changed: 114 additions & 159 deletions

‎src/gpu/compute_forces_acoustic_cuda.cu‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -137,9 +137,6 @@ TRACE("Kernel_2_acoustic");
137137
mp->d_wxgll,
138138
d_rhostore,
139139
mp->d_spec_to_pml,
140-
mp->ALPHA_MAX_PML,
141-
mp->d0_max_acoustic,
142-
mp->abscissa_norm,
143140
mp->nspec_pml_x,
144141
mp->nspec_pml_z,
145142
mp->deltat,

‎src/gpu/compute_forces_viscoelastic_cuda.cu‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -298,9 +298,6 @@ TRACE("Kernel_2");
298298
mp->simulation_type,
299299
mp->p_sv,
300300
mp->d_spec_to_pml,
301-
mp->ALPHA_MAX_PML,
302-
mp->d0_max_elastic,
303-
mp->abscissa_norm,
304301
mp->nspec_pml_x,
305302
mp->nspec_pml_z,
306303
mp->deltat,

‎src/gpu/kernels/Kernel_2_acoustic_impl.cu‎

Lines changed: 40 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -364,9 +364,6 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
364364
realw_const_p d_wxgll,
365365
const realw* d_rhostore,
366366
int* d_spec_to_pml,
367-
realw ALPHA_MAX_PML,
368-
realw d0,
369-
realw* abs_normalized,
370367
int NSPEC_PML_X,
371368
int NSPEC_PML_Z,
372369
realw deltat,
@@ -414,7 +411,7 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
414411
// PML
415412
int ispec_pml;
416413
int offset_pml,offset_local_pml;
417-
realw alpha1,beta1,alphax,betax,abs_norm;
414+
realw alpha1,beta1,alphax,betax,alphaz,betaz;
418415
realw c1,c2;
419416
realw r1,r2,r3,r4,r5,r6;
420417
realw rhol,kappal;
@@ -505,46 +502,53 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
505502
offset_local_pml = (ispec_pml-(NSPEC_PML_X + NSPEC_PML_Z))*NGLL2 + tx; // local pml elements in range [0,NSPEC_PML_XZ-1]
506503

507504
// coefficients
505+
alphax = alphax_store[offset_pml];
506+
betax = betax_store[offset_pml];
507+
alphaz = alphaz_store[offset_pml];
508+
betaz = betaz_store[offset_pml];
509+
508510
// note: see pml_init.F90, compare to routine define_PML_coefficients()
509511
// (starting around line 900);
510-
// assumes coefficients are calculated w/out PML_PARAMETER_ADJUSTMENT
511-
if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
512-
// in CPML_X_ONLY or in CPML_Z_ONLY region
513-
abs_norm = abs_normalized[offset_pml];
512+
// K_MIN_PML must be == 1.0 and K_MAX_PML == 1.0,
513+
// thus, K_x == K_z == K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER == 1
514+
if (ispec_pml < NSPEC_PML_X){
515+
// in CPML_X_ONLY
516+
//
514517
// for CPML_X_ONLY:
515-
// alpha1 == alpha_x == ALPHA_MAX * (1 - abscissa)
518+
// alpha1 == alpha_x
516519
// alpha_z == 0
517520
//
518521
// beta1 == beta_x == alpha_x + d_x / K_x
519522
// with d_x == d0_x / damping_change_factor_acoustic * abscissa_normalized**NPOWER
520523
// K_x == K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER
521524
// == alpha_x + (d0_x / damping_change_factor_acoustic * abscissa_normalized**NPOWER) /
522525
// (K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER)
523-
// note: damping_change_factor_acoustic must be == 0.5 for this implementation,
524-
// K_MIN_PML must be == 1.0 and K_MAX_PML == 1.0,
525-
// and NPOWER must be == 2 (see parameter defined in pml_init.F90)
526-
// Thus,
527-
// K_x == K_z == K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER == 1
528-
// d_x == d0_x / damping_change_factor_acoustic * abscissa_normalized**NPOWER == 2 * d0_x * abscissa**2
529-
// and it follows, that
530-
// beta1 == beta_x == alpha_x + (2 * d0_x * abscissa**2)
531526
// beta_z == 0
532527
//
528+
alpha1 = alphax;
529+
beta1 = betax;
530+
alphaz = 0.f;
531+
betaz = 0.f;
532+
} else if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
533+
// in CPML_Z_ONLY region
534+
//
533535
// for CPML_Z_ONLY:
534536
// alpha1 == alpha_z == ALPHA_MAX * (1 - abscissa)
535537
// alpha_x == 0
536538
//
537539
// beta1 == beta_z == alpha_z + d_z / K_z
538540
// == alpha_z + (2 * d0_z * abscissa**2)
539541
// beta_x == 0
540-
alpha1 = ALPHA_MAX_PML * (1.f - abs_norm);
541-
beta1 = alpha1 + 2.f * d0 * abs_norm * abs_norm; // instead of d0_x_left, d0_x_right, .. this takes d0_max
542+
alpha1 = alphaz;
543+
beta1 = betaz;
544+
alphax = 0.f;
545+
betax = 0.f;
542546
} else {
543547
// in CPML_XZ region
544-
alpha1 = alphaz_store[offset_local_pml];
545-
beta1 = betaz_store[offset_local_pml];
546-
alphax = alphax_store[offset_local_pml];
547-
betax = betax_store[offset_local_pml];
548+
alpha1 = alphaz;
549+
beta1 = betaz;
550+
//alphax = alphax;
551+
//betax = betax;
548552
}
549553

550554
// Update memory variables of derivatives
@@ -584,7 +588,7 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
584588
c2 = __expf(-0.5f * deltat * beta1);
585589

586590
coef0_1 = c1 * c1;
587-
if (abs(alpha1) > 0.00001){
591+
if (abs(alpha1) > 0.00001f){
588592
// coef1_zx_1 == (1 - c1)/alpha1
589593
// coef2_zx_1 == coef1 * c1
590594
coef1_1 = (1.f - c1) / alpha1;
@@ -597,7 +601,7 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
597601
}
598602

599603
coef0_2 = c2 * c2;
600-
if (abs(beta1) > 0.00001){
604+
if (abs(beta1) > 0.00001f){
601605
// coef1_zx_2 == (1 - c2)/beta1
602606
// coef2_zx_2 == coef1 * c2
603607
coef1_2 = (1.f - c2) / beta1;
@@ -615,7 +619,7 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
615619
realw c4 = __expf(-0.5f * deltat * alphax);
616620

617621
coef0_3 = c3 * c3;
618-
if (abs(betax) > 0.00001){
622+
if (abs(betax) > 0.00001f){
619623
// coef1_zx_2 == (1 - c3)/betax
620624
// coef2_zx_2 == coef1 * c3
621625
coef1_3 = (1.f - c3) / betax;
@@ -628,7 +632,7 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
628632
}
629633

630634
coef0_4 = c4 * c4;
631-
if (abs(alphax) > 0.00001){
635+
if (abs(alphax) > 0.00001f){
632636
// coef1_zx_1 == (1 - c4)/alphax
633637
// coef2_zx_1 == coef1 * c4
634638
coef1_4 = (1.f - c4) / alphax;
@@ -1006,8 +1010,15 @@ template __global__ void Kernel_2_acoustic_PML_impl<1>(const int,const int*,cons
10061010
realw_const_p,realw_p,
10071011
const realw*, const realw*,const realw*,const realw*,
10081012
realw_const_p,realw_const_p,realw_const_p,
1009-
const realw*,int*,realw,realw,realw*,
1010-
int,int,realw,
1013+
const realw*,int*,int,int,realw,
1014+
realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,
1015+
realw_p,realw*,realw*,realw*,realw*,realw*);
1016+
1017+
template __global__ void Kernel_2_acoustic_PML_impl<3>(const int,const int*,const int*,const int,const int,
1018+
realw_const_p,realw_p,
1019+
const realw*, const realw*,const realw*,const realw*,
1020+
realw_const_p,realw_const_p,realw_const_p,
1021+
const realw*,int*,int,int,realw,
10111022
realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,
10121023
realw_p,realw*,realw*,realw*,realw*,realw*);
10131024

‎src/gpu/kernels/Kernel_2_viscoelastic_impl.cu‎

Lines changed: 38 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -652,9 +652,6 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
652652
const int simulation_type,
653653
const int p_sv,
654654
const int* d_spec_to_pml,
655-
realw ALPHA_MAX_PML,
656-
realw d0,
657-
realw* abs_normalized,
658655
int NSPEC_PML_X,
659656
int NSPEC_PML_Z,
660657
realw deltat,
@@ -713,7 +710,7 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
713710
// PML
714711
int ispec_pml;
715712
int offset_pml,offset_local_pml;
716-
realw alpha1,beta1,alphax,betax,abs_norm;
713+
realw alpha1,beta1,alphax,betax,alphaz,betaz;
717714
realw c1,c2;
718715
realw r1,r2,r3,r4,r5,r6,r7,r8;
719716
realw r9_x,r9_z,r10_x,r10_z;
@@ -802,46 +799,53 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
802799
offset_local_pml = (ispec_pml-(NSPEC_PML_X + NSPEC_PML_Z))*NGLL2 + tx; // local pml elements in range [0,NSPEC_PML_XZ-1]
803800

804801
// coefficients
802+
alphax = alphax_store[offset_pml];
803+
betax = betax_store[offset_pml];
804+
alphaz = alphaz_store[offset_pml];
805+
betaz = betaz_store[offset_pml];
806+
805807
// note: see pml_init.F90, compare to routine define_PML_coefficients()
806808
// (starting around line 900);
807-
// assumes coefficients are calculated w/out PML_PARAMETER_ADJUSTMENT
808-
if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
809-
// in CPML_X_ONLY or in CPML_Z_ONLY region
810-
abs_norm = abs_normalized[offset_pml];
809+
// K_MIN_PML must be == 1.0 and K_MAX_PML == 1.0,
810+
// thus, K_x == K_z == K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER == 1
811+
if (ispec_pml < NSPEC_PML_X){
812+
// in CPML_X_ONLY
813+
//
811814
// for CPML_X_ONLY:
812-
// (see pml_init.F90, routine `define_PML_coefficients`)
813-
// alpha1 == alpha_x == ALPHA_MAX * (1 - abscissa)
815+
// alpha1 == alpha_x
814816
// alpha_z == 0
815817
//
816818
// beta1 == beta_x == alpha_x + d_x / K_x
817-
// == alpha_x + (d0_x / damping_change_factor_elastic * abscissa_normalized**NPOWER) /
819+
// with d_x == d0_x / damping_change_factor_acoustic * abscissa_normalized**NPOWER
820+
// K_x == K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER
821+
// == alpha_x + (d0_x / damping_change_factor_acoustic * abscissa_normalized**NPOWER) /
818822
// (K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER)
819-
// note: damping_change_factor_elastic must be == 1.0 for this implementation,
820-
// K_MIN_PML must be == 1.0 and K_MAX_PML == 1.0,
821-
// and NPOWER must be == 2 (see parameter defined in pml_init.F90)
822-
// Thus,
823-
// K_x == K_z == K_MIN_PML + (K_MAX_PML - 1.0d0) * abscissa_normalized**NPOWER == 1
824-
// d_x == d0_x / damping_change_factor_elastic * abscissa_normalized**NPOWER == d0_x * abscissa**2
825-
//
826-
// it follows, that
827-
// beta1 == beta_x == alpha_x + (d0_x * abscissa**2)
828823
// beta_z == 0
829824
//
825+
alpha1 = alphax;
826+
beta1 = betax;
827+
alphaz = 0.f;
828+
betaz = 0.f;
829+
} else if (ispec_pml < (NSPEC_PML_X + NSPEC_PML_Z)){
830+
// in CPML_Z_ONLY region
831+
//
830832
// for CPML_Z_ONLY:
831833
// alpha1 == alpha_z == ALPHA_MAX * (1 - abscissa)
832834
// alpha_x == 0
833835
//
834836
// beta1 == beta_z == alpha_z + d_z / K_z
835-
// == alpha_z + (d0_z * abscissa**2)
837+
// == alpha_z + (2 * d0_z * abscissa**2)
836838
// beta_x == 0
837-
alpha1 = ALPHA_MAX_PML * (1.f - abs_norm) ;
838-
beta1 = alpha1 + d0 * abs_norm * abs_norm; // instead of d0_x_left, d0_x_right, .. this takes d0_max
839-
} else{
839+
alpha1 = alphaz;
840+
beta1 = betaz;
841+
alphax = 0.f;
842+
betax = 0.f;
843+
} else {
840844
// in CPML_XZ region
841-
alpha1 = alphaz_store[offset_local_pml];
842-
beta1 = betaz_store[offset_local_pml];
843-
alphax = alphax_store[offset_local_pml];
844-
betax = betax_store[offset_local_pml];
845+
alpha1 = alphaz;
846+
beta1 = betaz;
847+
//alphax = alphax;
848+
//betax = betax;
845849
}
846850

847851
// Update memory variables of derivatives
@@ -880,7 +884,7 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
880884
c2 = __expf(-0.5f * deltat * beta1);
881885

882886
coef0_1 = c1 * c1;
883-
if (abs(alpha1) > 0.00001){
887+
if (abs(alpha1) > 0.00001f){
884888
// coef1_zx_1 == (1 - c1)/alpha1
885889
// coef2_zx_1 == coef1 * c1
886890
coef1_1 = (1.f - c1) / alpha1;
@@ -893,7 +897,7 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
893897
}
894898

895899
coef0_2 = c2 * c2;
896-
if (abs(beta1) > 0.00001){
900+
if (abs(beta1) > 0.00001f){
897901
// coef1_zx_2 == (1 - c2)/beta1
898902
// coef2_zx_2 == coef1 * c2
899903
coef1_2 = (1.f - c2) / beta1;
@@ -911,7 +915,7 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
911915
realw c4 = __expf(-0.5f * deltat * alphax);
912916

913917
coef0_3 = c3 * c3;
914-
if (abs(betax) > 0.00001){
918+
if (abs(betax) > 0.00001f){
915919
// coef1_zx_2 == (1 - c3)/betax
916920
// coef2_zx_2 == coef1 * c3
917921
coef1_3 = (1.f - c3) / betax;
@@ -924,7 +928,7 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
924928
}
925929

926930
coef0_4 = c4 * c4;
927-
if (abs(alphax) > 0.00001){
931+
if (abs(alphax) > 0.00001f){
928932
// coef1_zx_1 == (1 - c4)/alphax
929933
// coef2_zx_1 == coef1 * c4
930934
coef1_4 = (1.f - c4) / alphax;
@@ -2097,8 +2101,7 @@ template __global__ void Kernel_2_noatt_iso_PML_impl<1>(const int,const int*,con
20972101
realw_const_p,realw_p,
20982102
realw*,realw*,realw*,realw*,
20992103
realw_const_p,realw_const_p,realw_const_p,
2100-
realw*,realw*,const int,const int,const int*,
2101-
realw,realw,realw*,int,int,realw,
2104+
realw*,realw*,const int,const int,const int*,int,int,realw,
21022105
realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,
21032106
realw*,realw*,realw*,realw*,realw*,
21042107
realw_p,const realw*,
@@ -2108,8 +2111,7 @@ template __global__ void Kernel_2_noatt_iso_PML_impl<3>(const int,const int*,con
21082111
realw_const_p,realw_p,
21092112
realw*,realw*,realw*,realw*,
21102113
realw_const_p,realw_const_p,realw_const_p,
2111-
realw*,realw*,const int,const int,const int*,
2112-
realw,realw,realw*,int,int,realw,
2114+
realw*,realw*,const int,const int,const int*,int,int,realw,
21132115
realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,realw*,
21142116
realw*,realw*,realw*,realw*,realw*,
21152117
realw_p,const realw*,

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

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -454,9 +454,6 @@ Kernel_2_acoustic_PML_impl(const int nb_blocks_to_compute,
454454
realw_const_p d_wxgll,
455455
const realw* d_rhostore,
456456
int* d_spec_to_pml,
457-
realw ALPHA_MAX_PML,
458-
realw d0,
459-
realw* abs_normalized,
460457
int NSPEC_PML_X,
461458
int NSPEC_PML_Z,
462459
realw deltat,
@@ -537,9 +534,6 @@ Kernel_2_noatt_iso_PML_impl(const int nb_blocks_to_compute,
537534
const int simulation_type,
538535
const int p_sv,
539536
const int* d_spec_to_pml,
540-
realw ALPHA_MAX_PML,
541-
realw d0,
542-
realw* abs_normalized,
543537
int NSPEC_PML_X,
544538
int NSPEC_PML_Z,
545539
realw deltat,

‎src/gpu/mesh_constants_cuda.h‎

Lines changed: 0 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -575,12 +575,7 @@ typedef struct mesh_ {
575575
realw* betax_store;
576576
realw* betaz_store;
577577

578-
realw* abscissa_norm;
579-
realw ALPHA_MAX_PML;
580-
581578
// acoustic
582-
realw d0_max_acoustic;
583-
584579
realw* PML_dpotentialdxl_old;
585580
realw* PML_dpotentialdzl_old;
586581
realw* d_potential_old;
@@ -596,8 +591,6 @@ typedef struct mesh_ {
596591
int* d_pml_abs_points_acoustic;
597592

598593
// elastic
599-
realw d0_max_elastic;
600-
601594
realw* PML_dux_dxl_old;
602595
realw* PML_dux_dzl_old;
603596
realw* PML_duz_dxl_old;

0 commit comments

Comments
 (0)