[cig-commits] r20532 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: UTILS src/cuda src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Jul 20 14:52:01 PDT 2012
Author: danielpeter
Date: 2012-07-20 14:52:01 -0700 (Fri, 20 Jul 2012)
New Revision: 20532
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/UTILS/update_headers_change_word_f90.pl
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_kernels_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
Log:
updates cuda routines to account for 1D and 3D attenuation arrays
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/UTILS/update_headers_change_word_f90.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/UTILS/update_headers_change_word_f90.pl 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/UTILS/update_headers_change_word_f90.pl 2012-07-20 21:52:01 UTC (rev 20532)
@@ -42,6 +42,9 @@
# suppress trailing white spaces and carriage return
$line =~ s/\s*$//;
+# converts tabs to space
+ $line =~ s/\t/ /g;
+
# change the version number and copyright information
# $line =~ s#\(c\) California Institute of Technology and University of Pau, October 2007#\(c\) California Institute of Technology and University of Pau, November 2007#og;
# $line =~ s#rmass_sigma#rmass_time_integral_of_sigma#og;
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu 2012-07-20 21:52:01 UTC (rev 20532)
@@ -313,7 +313,8 @@
reald epsilondev_xx_loc,reald epsilondev_yy_loc,reald epsilondev_xy_loc,
reald epsilondev_xz_loc,reald epsilondev_yz_loc,
int ANISOTROPY,
- realw* d_c44store
+ realw* d_c44store,
+ int ATTENUATION_3D
){
int i_sls;
@@ -346,8 +347,12 @@
// either mustore(i,j,k,ispec) * factor_common(i_sls,i,j,k,ispec)
// or factor_common(i_sls,:,:,:,ispec) * c44store(:,:,:,ispec)
- factor_loc = fac * factor_common[offset];
-
+ if( ATTENUATION_3D ){
+ factor_loc = fac * factor_common[offset];
+ }else{
+ factor_loc = fac * factor_common[i_sls + N_SLS*working_element];
+ }
+
alphaval_loc = alphaval[i_sls]; // (i_sls)
betaval_loc = betaval[i_sls];
gammaval_loc = gammaval[i_sls];
@@ -917,6 +922,7 @@
int ATTENUATION,
int ATTENUATION_NEW,
int USE_ATTENUATION_MIMIC,
+ int ATTENUATION_3D,
realw* one_minus_sum_beta,realw* factor_common,
realw* R_xx, realw* R_yy, realw* R_xy, realw* R_xz, realw* R_yz,
realw* alphaval,realw* betaval,realw* gammaval,
@@ -1034,26 +1040,26 @@
#endif
if(ATTENUATION){
- if(ATTENUATION_NEW){
- // takes new routines
- // use first order Taylor expansion of displacement for local storage of stresses
- // at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION_NEW){
+ // takes new routines
+ // use first order Taylor expansion of displacement for local storage of stresses
+ // at this current time step, to fix attenuation in a consistent way
#ifdef USE_TEXTURES
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
#else
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
#endif
- }
- else{
- // takes old routines
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx];
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx];
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx];
- }
+ }
+ else{
+ // takes old routines
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx];
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx];
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx];
+ }
}
}
@@ -1103,40 +1109,39 @@
if( ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
+ tempx1l_att = 0.f;
+ tempx2l_att = 0.f;
+ tempx3l_att = 0.f;
- tempx1l_att = 0.f;
- tempx2l_att = 0.f;
- tempx3l_att = 0.f;
+ tempy1l_att = 0.f;
+ tempy2l_att = 0.f;
+ tempy3l_att = 0.f;
- tempy1l_att = 0.f;
- tempy2l_att = 0.f;
- tempy3l_att = 0.f;
+ tempz1l_att = 0.f;
+ tempz2l_att = 0.f;
+ tempz3l_att = 0.f;
- tempz1l_att = 0.f;
- tempz2l_att = 0.f;
- tempz3l_att = 0.f;
+ for (l=0;l<NGLLX;l++) {
+ hp1 = d_hprime_xx[l*NGLLX+I];
+ offset = K*NGLL2+J*NGLLX+l;
+ tempx1l_att += s_dummyx_loc_att[offset]*hp1;
+ tempy1l_att += s_dummyy_loc_att[offset]*hp1;
+ tempz1l_att += s_dummyz_loc_att[offset]*hp1;
- for (l=0;l<NGLLX;l++) {
- hp1 = d_hprime_xx[l*NGLLX+I];
- offset = K*NGLL2+J*NGLLX+l;
- tempx1l_att += s_dummyx_loc_att[offset]*hp1;
- tempy1l_att += s_dummyy_loc_att[offset]*hp1;
- tempz1l_att += s_dummyz_loc_att[offset]*hp1;
+ hp2 = d_hprime_xx[l*NGLLX+J];
+ offset = K*NGLL2+l*NGLLX+I;
+ tempx2l_att += s_dummyx_loc_att[offset]*hp2;
+ tempy2l_att += s_dummyy_loc_att[offset]*hp2;
+ tempz2l_att += s_dummyz_loc_att[offset]*hp2;
- hp2 = d_hprime_xx[l*NGLLX+J];
- offset = K*NGLL2+l*NGLLX+I;
- tempx2l_att += s_dummyx_loc_att[offset]*hp2;
- tempy2l_att += s_dummyy_loc_att[offset]*hp2;
- tempz2l_att += s_dummyz_loc_att[offset]*hp2;
+ hp3 = d_hprime_xx[l*NGLLX+K];
+ offset = l*NGLL2+J*NGLLX+I;
+ tempx3l_att += s_dummyx_loc_att[offset]*hp3;
+ tempy3l_att += s_dummyy_loc_att[offset]*hp3;
+ tempz3l_att += s_dummyz_loc_att[offset]*hp3;
- hp3 = d_hprime_xx[l*NGLLX+K];
- offset = l*NGLL2+J*NGLLX+I;
- tempx3l_att += s_dummyx_loc_att[offset]*hp3;
- tempy3l_att += s_dummyy_loc_att[offset]*hp3;
- tempz3l_att += s_dummyz_loc_att[offset]*hp3;
-
- }
+ }
}
#else
@@ -1195,61 +1200,60 @@
+ s_dummyz_loc[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
if( ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
+ tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
-
- tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
}
#endif
@@ -1288,66 +1292,69 @@
duzdyl_plus_duydzl = duzdyl + duydzl;
if( ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
+ duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
+ duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
- duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
- duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
- duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
+ duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
+ duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
+ duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
- duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
- duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
- duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
+ duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
+ duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
+ duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
- duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
- duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
- duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
+ // precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
- // precompute some sums to save CPU time
- duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
- duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
- duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl_att - templ;
+ epsilondev_yy_loc = duydyl_att - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl_att - templ;
- epsilondev_yy_loc = duydyl_att - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
-
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
- epsilon_trace_over_3[tx] = templ;
- }else{
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
+ epsilon_trace_over_3[tx] = templ;
+ }else{
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}else{
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl - templ;
- epsilondev_yy_loc = duydyl - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl - templ;
+ epsilondev_yy_loc = duydyl - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
- epsilon_trace_over_3[tx] = templ;
- }else{
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
+ epsilon_trace_over_3[tx] = templ;
+ }else{
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}
// attenuation
if(ATTENUATION){
// use unrelaxed parameters if attenuation
- one_minus_sum_beta_use = one_minus_sum_beta[tx+working_element*NGLL3]; // (i,j,k,ispec)
+ if( ATTENUATION_3D ){
+ one_minus_sum_beta_use = one_minus_sum_beta[tx+working_element*NGLL3]; // (i,j,k,ispec)
+ }else{
+ one_minus_sum_beta_use = one_minus_sum_beta[working_element]; // (1,1,1,ispec)
+ }
minus_sum_beta = one_minus_sum_beta_use - 1.0f;
}
@@ -1603,7 +1610,7 @@
R_xx,R_yy,R_xy,R_xz,R_yz,
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz,
epsilondev_xx_loc,epsilondev_yy_loc,epsilondev_xy_loc,epsilondev_xz_loc,epsilondev_yz_loc,
- ANISOTROPY,d_c44store);
+ ANISOTROPY,d_c44store,ATTENUATION_3D);
}
// save deviatoric strain for Runge-Kutta scheme
@@ -1725,6 +1732,7 @@
mp->attenuation,
mp->attenuation_new,
mp->use_attenuation_mimic,
+ mp->attenuation_3D,
d_one_minus_sum_beta,d_factor_common,
d_R_xx,d_R_yy,d_R_xy,d_R_xz,d_R_yz,
mp->d_alphaval,mp->d_betaval,mp->d_gammaval,
@@ -1772,6 +1780,7 @@
mp->attenuation,
mp->attenuation_new,
mp->use_attenuation_mimic,
+ mp->attenuation_3D,
d_one_minus_sum_beta,d_factor_common,
d_b_R_xx,d_b_R_yy,d_b_R_xy,d_b_R_xz,d_b_R_yz,
mp->d_b_alphaval,mp->d_b_betaval,mp->d_b_gammaval,
@@ -1864,7 +1873,11 @@
// array offsets
color_offset = (mp->nspec_outer_crust_mantle) * NGLL3_PADDED;
color_offset_nonpadded = (mp->nspec_outer_crust_mantle) * NGLL3;
- color_offset_nonpadded_att2 = (mp->nspec_outer_crust_mantle) * NGLL3 * N_SLS;
+ if( mp->attenuation_3D ){
+ color_offset_nonpadded_att2 = (mp->nspec_outer_crust_mantle) * NGLL3 * N_SLS;
+ }else{
+ color_offset_nonpadded_att2 = (mp->nspec_outer_crust_mantle) * 1 * N_SLS;
+ }
color_offset_ispec = mp->nspec_outer_crust_mantle;
}
@@ -1950,7 +1963,11 @@
// for no-aligned arrays
color_offset_nonpadded += nb_blocks_to_compute * NGLL3;
// for factor_common array
- color_offset_nonpadded_att2 += nb_blocks_to_compute * NGLL3 * N_SLS;
+ if( mp->attenuation_3D ){
+ color_offset_nonpadded_att2 += nb_blocks_to_compute * NGLL3 * N_SLS;
+ }else{
+ color_offset_nonpadded_att2 += nb_blocks_to_compute * 1 * N_SLS;
+ }
// for array(ispec)
color_offset_ispec += nb_blocks_to_compute;
}
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu 2012-07-20 21:52:01 UTC (rev 20532)
@@ -92,7 +92,8 @@
realw* epsilondev_xx,realw* epsilondev_yy,realw* epsilondev_xy,
realw* epsilondev_xz,realw* epsilondev_yz,
reald epsilondev_xx_loc,reald epsilondev_yy_loc,reald epsilondev_xy_loc,
- reald epsilondev_xz_loc,reald epsilondev_yz_loc
+ reald epsilondev_xz_loc,reald epsilondev_yz_loc,
+ int ATTENUATION_3D
){
int i_sls;
@@ -117,8 +118,11 @@
// index for (i_sls,i,j,k,ispec)
offset = i_sls + N_SLS*(tx + NGLL3*working_element);
- factor_loc = mul * factor_common[offset]; //mustore(i,j,k,ispec) * factor_common(i_sls,i,j,k,ispec)
-
+ if( ATTENUATION_3D ){
+ factor_loc = mul * factor_common[offset]; //mustore(i,j,k,ispec) * factor_common(i_sls,i,j,k,ispec)
+ }else{
+ factor_loc = mul * factor_common[i_sls + N_SLS*working_element]; //mustore(i,j,k,ispec) * factor_common(i_sls,1,1,1,ispec)
+ }
alphaval_loc = alphaval[i_sls]; // (i_sls)
betaval_loc = betaval[i_sls];
gammaval_loc = gammaval[i_sls];
@@ -292,45 +296,46 @@
/* ----------------------------------------------------------------------------------------------- */
__global__ void Kernel_2_inner_core_impl(int nb_blocks_to_compute,
- int NGLOB,
- int* d_ibool,
- int* d_idoubling,
- int* d_phase_ispec_inner,
- int num_phase_ispec,
- int d_iphase,
- realw d_deltat,
- int use_mesh_coloring_gpu,
- realw* d_displ,
- realw* d_veloc,
- realw* d_accel,
- realw* d_xix, realw* d_xiy, realw* d_xiz,
- realw* d_etax, realw* d_etay, realw* d_etaz,
- realw* d_gammax, realw* d_gammay, realw* d_gammaz,
- realw* d_hprime_xx, realw* d_hprime_yy, realw* d_hprime_zz,
- realw* d_hprimewgll_xx, realw* d_hprimewgll_yy, realw* d_hprimewgll_zz,
- realw* d_wgllwgll_xy,realw* d_wgllwgll_xz,realw* d_wgllwgll_yz,
- realw* d_kappav,
- realw* d_muv,
- int COMPUTE_AND_STORE_STRAIN,
- realw* epsilondev_xx,realw* epsilondev_yy,realw* epsilondev_xy,
- realw* epsilondev_xz,realw* epsilondev_yz,
- realw* epsilon_trace_over_3,
- int SIMULATION_TYPE,
- int ATTENUATION,
- int ATTENUATION_NEW,
- int USE_ATTENUATION_MIMIC,
- realw* one_minus_sum_beta,realw* factor_common,
- realw* R_xx, realw* R_yy, realw* R_xy, realw* R_xz, realw* R_yz,
- realw* alphaval,realw* betaval,realw* gammaval,
- int ANISOTROPY,
- realw* d_c11store,realw* d_c12store,realw* d_c13store,
- realw* d_c33store,realw* d_c44store,
- int GRAVITY,
- realw* d_xstore,realw* d_ystore,realw* d_zstore,
- realw* d_minus_gravity_table,
- realw* d_minus_deriv_gravity_table,
- realw* d_density_table,
- realw* wgll_cube){
+ int NGLOB,
+ int* d_ibool,
+ int* d_idoubling,
+ int* d_phase_ispec_inner,
+ int num_phase_ispec,
+ int d_iphase,
+ realw d_deltat,
+ int use_mesh_coloring_gpu,
+ realw* d_displ,
+ realw* d_veloc,
+ realw* d_accel,
+ realw* d_xix, realw* d_xiy, realw* d_xiz,
+ realw* d_etax, realw* d_etay, realw* d_etaz,
+ realw* d_gammax, realw* d_gammay, realw* d_gammaz,
+ realw* d_hprime_xx, realw* d_hprime_yy, realw* d_hprime_zz,
+ realw* d_hprimewgll_xx, realw* d_hprimewgll_yy, realw* d_hprimewgll_zz,
+ realw* d_wgllwgll_xy,realw* d_wgllwgll_xz,realw* d_wgllwgll_yz,
+ realw* d_kappav,
+ realw* d_muv,
+ int COMPUTE_AND_STORE_STRAIN,
+ realw* epsilondev_xx,realw* epsilondev_yy,realw* epsilondev_xy,
+ realw* epsilondev_xz,realw* epsilondev_yz,
+ realw* epsilon_trace_over_3,
+ int SIMULATION_TYPE,
+ int ATTENUATION,
+ int ATTENUATION_NEW,
+ int USE_ATTENUATION_MIMIC,
+ int ATTENUATION_3D,
+ realw* one_minus_sum_beta,realw* factor_common,
+ realw* R_xx, realw* R_yy, realw* R_xy, realw* R_xz, realw* R_yz,
+ realw* alphaval,realw* betaval,realw* gammaval,
+ int ANISOTROPY,
+ realw* d_c11store,realw* d_c12store,realw* d_c13store,
+ realw* d_c33store,realw* d_c44store,
+ int GRAVITY,
+ realw* d_xstore,realw* d_ystore,realw* d_zstore,
+ realw* d_minus_gravity_table,
+ realw* d_minus_deriv_gravity_table,
+ realw* d_density_table,
+ realw* wgll_cube){
/* int bx = blockIdx.y*blockDim.x+blockIdx.x; //possible bug in original code*/
int bx = blockIdx.y*gridDim.x+blockIdx.x;
@@ -434,28 +439,28 @@
s_dummyz_loc[tx] = d_displ[iglob*3 + 2];
#endif
- if(ATTENUATION){
- if(ATTENUATION_NEW){
- // takes new routines
- // use first order Taylor expansion of displacement for local storage of stresses
- // at this current time step, to fix attenuation in a consistent way
+ if(ATTENUATION){
+ if(ATTENUATION_NEW){
+ // takes new routines
+ // use first order Taylor expansion of displacement for local storage of stresses
+ // at this current time step, to fix attenuation in a consistent way
#ifdef USE_TEXTURES
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
#else
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
#endif
- }
- else{
- // takes old routines
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx];
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx];
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx];
- }
- }
+ }
+ else{
+ // takes old routines
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx];
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx];
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx];
+ }
+ }
}
}
@@ -504,40 +509,39 @@
}
if( ATTENUATION ){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
+ tempx1l_att = 0.f;
+ tempx2l_att = 0.f;
+ tempx3l_att = 0.f;
- tempx1l_att = 0.f;
- tempx2l_att = 0.f;
- tempx3l_att = 0.f;
+ tempy1l_att = 0.f;
+ tempy2l_att = 0.f;
+ tempy3l_att = 0.f;
- tempy1l_att = 0.f;
- tempy2l_att = 0.f;
- tempy3l_att = 0.f;
+ tempz1l_att = 0.f;
+ tempz2l_att = 0.f;
+ tempz3l_att = 0.f;
- tempz1l_att = 0.f;
- tempz2l_att = 0.f;
- tempz3l_att = 0.f;
+ for (l=0;l<NGLLX;l++) {
+ hp1 = d_hprime_xx[l*NGLLX+I];
+ offset = K*NGLL2+J*NGLLX+l;
+ tempx1l_att += s_dummyx_loc_att[offset]*hp1;
+ tempy1l_att += s_dummyy_loc_att[offset]*hp1;
+ tempz1l_att += s_dummyz_loc_att[offset]*hp1;
- for (l=0;l<NGLLX;l++) {
- hp1 = d_hprime_xx[l*NGLLX+I];
- offset = K*NGLL2+J*NGLLX+l;
- tempx1l_att += s_dummyx_loc_att[offset]*hp1;
- tempy1l_att += s_dummyy_loc_att[offset]*hp1;
- tempz1l_att += s_dummyz_loc_att[offset]*hp1;
+ hp2 = d_hprime_xx[l*NGLLX+J];
+ offset = K*NGLL2+l*NGLLX+I;
+ tempx2l_att += s_dummyx_loc_att[offset]*hp2;
+ tempy2l_att += s_dummyy_loc_att[offset]*hp2;
+ tempz2l_att += s_dummyz_loc_att[offset]*hp2;
- hp2 = d_hprime_xx[l*NGLLX+J];
- offset = K*NGLL2+l*NGLLX+I;
- tempx2l_att += s_dummyx_loc_att[offset]*hp2;
- tempy2l_att += s_dummyy_loc_att[offset]*hp2;
- tempz2l_att += s_dummyz_loc_att[offset]*hp2;
+ hp3 = d_hprime_xx[l*NGLLX+K];
+ offset = l*NGLL2+J*NGLLX+I;
+ tempx3l_att += s_dummyx_loc_att[offset]*hp3;
+ tempy3l_att += s_dummyy_loc_att[offset]*hp3;
+ tempz3l_att += s_dummyz_loc_att[offset]*hp3;
- hp3 = d_hprime_xx[l*NGLLX+K];
- offset = l*NGLL2+J*NGLLX+I;
- tempx3l_att += s_dummyx_loc_att[offset]*hp3;
- tempy3l_att += s_dummyy_loc_att[offset]*hp3;
- tempz3l_att += s_dummyz_loc_att[offset]*hp3;
-
- }
+ }
}
#else
@@ -596,61 +600,60 @@
+ s_dummyz_loc[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
if( ATTENUATION ){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
+ tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
-
- tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
}
#endif
@@ -689,56 +692,55 @@
duzdyl_plus_duydzl = duzdyl + duydzl;
if(ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
+ duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
+ duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
- duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
- duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
- duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
+ duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
+ duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
+ duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
- duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
- duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
- duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
+ duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
+ duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
+ duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
- duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
- duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
- duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
+ // precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
- // precompute some sums to save CPU time
- duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
- duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
- duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl_att - templ;
+ epsilondev_yy_loc = duydyl_att - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl_att - templ;
- epsilondev_yy_loc = duydyl_att - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
-
- if(SIMULATION_TYPE == 3) {
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(SIMULATION_TYPE == 3) {
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}else{
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl - templ;
- epsilondev_yy_loc = duydyl - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl - templ;
+ epsilondev_yy_loc = duydyl - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
- if(SIMULATION_TYPE == 3) {
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(SIMULATION_TYPE == 3) {
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}
// compute elements with an elastic isotropic rheology
@@ -748,8 +750,13 @@
// attenuation
if(ATTENUATION){
// use unrelaxed parameters if attenuation
- mul_iso = mul * one_minus_sum_beta[tx+working_element*NGLL3]; // (i,j,k,ispec)
- mul_aniso = mul *( one_minus_sum_beta[tx+working_element*NGLL3] - 1.0f );
+ if( ATTENUATION_3D ){
+ mul_iso = mul * one_minus_sum_beta[tx+working_element*NGLL3]; // (i,j,k,ispec)
+ mul_aniso = mul *( one_minus_sum_beta[tx+working_element*NGLL3] - 1.0f );
+ }else{
+ mul_iso = mul * one_minus_sum_beta[working_element]; // (1,1,1,ispec)
+ mul_aniso = mul *( one_minus_sum_beta[working_element] - 1.0f );
+ }
}else{
mul_iso = mul;
}
@@ -1017,7 +1024,8 @@
factor_common,alphaval,betaval,gammaval,
R_xx,R_yy,R_xy,R_xz,R_yz,
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz,
- epsilondev_xx_loc,epsilondev_yy_loc,epsilondev_xy_loc,epsilondev_xz_loc,epsilondev_yz_loc);
+ epsilondev_xx_loc,epsilondev_yy_loc,epsilondev_xy_loc,epsilondev_xz_loc,epsilondev_yz_loc,
+ ATTENUATION_3D);
}
// save deviatoric strain for Runge-Kutta scheme
@@ -1044,42 +1052,42 @@
/* ----------------------------------------------------------------------------------------------- */
void Kernel_2_inner_core(int nb_blocks_to_compute,Mesh* mp,
- realw d_deltat,
- int d_iphase,
- int* d_ibool,
- int* d_idoubling,
- realw* d_xix,realw* d_xiy,realw* d_xiz,
- realw* d_etax,realw* d_etay,realw* d_etaz,
- realw* d_gammax,realw* d_gammay,realw* d_gammaz,
- realw* d_kappav,
- realw* d_muv,
- realw* d_epsilondev_xx,
- realw* d_epsilondev_yy,
- realw* d_epsilondev_xy,
- realw* d_epsilondev_xz,
- realw* d_epsilondev_yz,
- realw* d_epsilon_trace_over_3,
- realw* d_one_minus_sum_beta,
- realw* d_factor_common,
- realw* d_R_xx,
- realw* d_R_yy,
- realw* d_R_xy,
- realw* d_R_xz,
- realw* d_R_yz,
- realw* d_b_epsilondev_xx,
- realw* d_b_epsilondev_yy,
- realw* d_b_epsilondev_xy,
- realw* d_b_epsilondev_xz,
- realw* d_b_epsilondev_yz,
- realw* d_b_epsilon_trace_over_3,
- realw* d_b_R_xx,
- realw* d_b_R_yy,
- realw* d_b_R_xy,
- realw* d_b_R_xz,
- realw* d_b_R_yz,
- realw* d_c11store,realw* d_c12store,realw* d_c13store,
- realw* d_c33store,realw* d_c44store
- ){
+ realw d_deltat,
+ int d_iphase,
+ int* d_ibool,
+ int* d_idoubling,
+ realw* d_xix,realw* d_xiy,realw* d_xiz,
+ realw* d_etax,realw* d_etay,realw* d_etaz,
+ realw* d_gammax,realw* d_gammay,realw* d_gammaz,
+ realw* d_kappav,
+ realw* d_muv,
+ realw* d_epsilondev_xx,
+ realw* d_epsilondev_yy,
+ realw* d_epsilondev_xy,
+ realw* d_epsilondev_xz,
+ realw* d_epsilondev_yz,
+ realw* d_epsilon_trace_over_3,
+ realw* d_one_minus_sum_beta,
+ realw* d_factor_common,
+ realw* d_R_xx,
+ realw* d_R_yy,
+ realw* d_R_xy,
+ realw* d_R_xz,
+ realw* d_R_yz,
+ realw* d_b_epsilondev_xx,
+ realw* d_b_epsilondev_yy,
+ realw* d_b_epsilondev_xy,
+ realw* d_b_epsilondev_xz,
+ realw* d_b_epsilondev_yz,
+ realw* d_b_epsilon_trace_over_3,
+ realw* d_b_R_xx,
+ realw* d_b_R_yy,
+ realw* d_b_R_xy,
+ realw* d_b_R_xz,
+ realw* d_b_R_yz,
+ realw* d_c11store,realw* d_c12store,realw* d_c13store,
+ realw* d_c33store,realw* d_c44store
+ ){
#ifdef ENABLE_VERY_SLOW_ERROR_CHECKING
exit_on_cuda_error("before kernel Kernel_2_inner_core");
@@ -1108,94 +1116,96 @@
// cudaEventRecord( start, 0 );
Kernel_2_inner_core_impl<<<grid,threads>>>(nb_blocks_to_compute,
- mp->NGLOB_INNER_CORE,
- d_ibool,
- d_idoubling,
- mp->d_phase_ispec_inner_inner_core,
- mp->num_phase_ispec_inner_core,
- d_iphase,
- d_deltat,
- mp->use_mesh_coloring_gpu,
- mp->d_displ_inner_core,
- mp->d_veloc_inner_core,
- mp->d_accel_inner_core,
- d_xix, d_xiy, d_xiz,
- d_etax, d_etay, d_etaz,
- d_gammax, d_gammay, d_gammaz,
- mp->d_hprime_xx, mp->d_hprime_yy, mp->d_hprime_zz,
- mp->d_hprimewgll_xx, mp->d_hprimewgll_yy, mp->d_hprimewgll_zz,
- mp->d_wgllwgll_xy, mp->d_wgllwgll_xz, mp->d_wgllwgll_yz,
- d_kappav, d_muv,
- mp->compute_and_store_strain,
- d_epsilondev_xx,
- d_epsilondev_yy,
- d_epsilondev_xy,
- d_epsilondev_xz,
- d_epsilondev_yz,
- d_epsilon_trace_over_3,
- mp->simulation_type,
- mp->attenuation,
- mp->attenuation_new,
- mp->use_attenuation_mimic,
- d_one_minus_sum_beta,
- d_factor_common,
- d_R_xx,d_R_yy,d_R_xy,d_R_xz,d_R_yz,
- mp->d_alphaval,mp->d_betaval,mp->d_gammaval,
- mp->anisotropic_inner_core,
- d_c11store,d_c12store,d_c13store,
- d_c33store,d_c44store,
- mp->gravity,
- mp->d_xstore_inner_core,mp->d_ystore_inner_core,mp->d_zstore_inner_core,
- mp->d_minus_gravity_table,
- mp->d_minus_deriv_gravity_table,
- mp->d_density_table,
- mp->d_wgll_cube);
+ mp->NGLOB_INNER_CORE,
+ d_ibool,
+ d_idoubling,
+ mp->d_phase_ispec_inner_inner_core,
+ mp->num_phase_ispec_inner_core,
+ d_iphase,
+ d_deltat,
+ mp->use_mesh_coloring_gpu,
+ mp->d_displ_inner_core,
+ mp->d_veloc_inner_core,
+ mp->d_accel_inner_core,
+ d_xix, d_xiy, d_xiz,
+ d_etax, d_etay, d_etaz,
+ d_gammax, d_gammay, d_gammaz,
+ mp->d_hprime_xx, mp->d_hprime_yy, mp->d_hprime_zz,
+ mp->d_hprimewgll_xx, mp->d_hprimewgll_yy, mp->d_hprimewgll_zz,
+ mp->d_wgllwgll_xy, mp->d_wgllwgll_xz, mp->d_wgllwgll_yz,
+ d_kappav, d_muv,
+ mp->compute_and_store_strain,
+ d_epsilondev_xx,
+ d_epsilondev_yy,
+ d_epsilondev_xy,
+ d_epsilondev_xz,
+ d_epsilondev_yz,
+ d_epsilon_trace_over_3,
+ mp->simulation_type,
+ mp->attenuation,
+ mp->attenuation_new,
+ mp->use_attenuation_mimic,
+ mp->attenuation_3D,
+ d_one_minus_sum_beta,
+ d_factor_common,
+ d_R_xx,d_R_yy,d_R_xy,d_R_xz,d_R_yz,
+ mp->d_alphaval,mp->d_betaval,mp->d_gammaval,
+ mp->anisotropic_inner_core,
+ d_c11store,d_c12store,d_c13store,
+ d_c33store,d_c44store,
+ mp->gravity,
+ mp->d_xstore_inner_core,mp->d_ystore_inner_core,mp->d_zstore_inner_core,
+ mp->d_minus_gravity_table,
+ mp->d_minus_deriv_gravity_table,
+ mp->d_density_table,
+ mp->d_wgll_cube);
if(mp->simulation_type == 3) {
Kernel_2_inner_core_impl<<< grid,threads>>>(nb_blocks_to_compute,
- mp->NGLOB_INNER_CORE,
- d_ibool,
- d_idoubling,
- mp->d_phase_ispec_inner_inner_core,
- mp->num_phase_ispec_inner_core,
- d_iphase,
- d_deltat,
- mp->use_mesh_coloring_gpu,
- mp->d_b_displ_inner_core,
- mp->d_b_veloc_inner_core,
- mp->d_b_accel_inner_core,
- d_xix, d_xiy, d_xiz,
- d_etax, d_etay, d_etaz,
- d_gammax, d_gammay, d_gammaz,
- mp->d_hprime_xx, mp->d_hprime_yy, mp->d_hprime_zz,
- mp->d_hprimewgll_xx, mp->d_hprimewgll_yy, mp->d_hprimewgll_zz,
- mp->d_wgllwgll_xy, mp->d_wgllwgll_xz, mp->d_wgllwgll_yz,
- d_kappav, d_muv,
- mp->compute_and_store_strain,
- d_b_epsilondev_xx,
- d_b_epsilondev_yy,
- d_b_epsilondev_xy,
- d_b_epsilondev_xz,
- d_b_epsilondev_yz,
- d_b_epsilon_trace_over_3,
- mp->simulation_type,
- mp->attenuation,
- mp->attenuation_new,
- mp->use_attenuation_mimic,
- d_one_minus_sum_beta,
- d_factor_common,
- d_b_R_xx,d_b_R_yy,d_b_R_xy,d_b_R_xz,d_b_R_yz,
- mp->d_b_alphaval,mp->d_b_betaval,mp->d_b_gammaval,
- mp->anisotropic_inner_core,
- d_c11store,d_c12store,d_c13store,
- d_c33store,d_c44store,
- mp->gravity,
- mp->d_xstore_inner_core,mp->d_ystore_inner_core,mp->d_zstore_inner_core,
- mp->d_minus_gravity_table,
- mp->d_minus_deriv_gravity_table,
- mp->d_density_table,
- mp->d_wgll_cube);
+ mp->NGLOB_INNER_CORE,
+ d_ibool,
+ d_idoubling,
+ mp->d_phase_ispec_inner_inner_core,
+ mp->num_phase_ispec_inner_core,
+ d_iphase,
+ d_deltat,
+ mp->use_mesh_coloring_gpu,
+ mp->d_b_displ_inner_core,
+ mp->d_b_veloc_inner_core,
+ mp->d_b_accel_inner_core,
+ d_xix, d_xiy, d_xiz,
+ d_etax, d_etay, d_etaz,
+ d_gammax, d_gammay, d_gammaz,
+ mp->d_hprime_xx, mp->d_hprime_yy, mp->d_hprime_zz,
+ mp->d_hprimewgll_xx, mp->d_hprimewgll_yy, mp->d_hprimewgll_zz,
+ mp->d_wgllwgll_xy, mp->d_wgllwgll_xz, mp->d_wgllwgll_yz,
+ d_kappav, d_muv,
+ mp->compute_and_store_strain,
+ d_b_epsilondev_xx,
+ d_b_epsilondev_yy,
+ d_b_epsilondev_xy,
+ d_b_epsilondev_xz,
+ d_b_epsilondev_yz,
+ d_b_epsilon_trace_over_3,
+ mp->simulation_type,
+ mp->attenuation,
+ mp->attenuation_new,
+ mp->use_attenuation_mimic,
+ mp->attenuation_3D,
+ d_one_minus_sum_beta,
+ d_factor_common,
+ d_b_R_xx,d_b_R_yy,d_b_R_xy,d_b_R_xz,d_b_R_yz,
+ mp->d_b_alphaval,mp->d_b_betaval,mp->d_b_gammaval,
+ mp->anisotropic_inner_core,
+ d_c11store,d_c12store,d_c13store,
+ d_c33store,d_c44store,
+ mp->gravity,
+ mp->d_xstore_inner_core,mp->d_ystore_inner_core,mp->d_zstore_inner_core,
+ mp->d_minus_gravity_table,
+ mp->d_minus_deriv_gravity_table,
+ mp->d_density_table,
+ mp->d_wgll_cube);
}
// cudaEventRecord( stop, 0 );
@@ -1218,8 +1228,8 @@
extern "C"
void FC_FUNC_(compute_forces_inner_core_cuda,
COMPUTE_FORCES_INNER_CORE_CUDA)(long* Mesh_pointer_f,
- realw* deltat,
- int* iphase) {
+ realw* deltat,
+ int* iphase) {
TRACE("compute_forces_inner_core_cuda");
@@ -1270,7 +1280,11 @@
// array offsets
color_offset = (mp->nspec_outer_inner_core) * NGLL3_PADDED;
color_offset_nonpadded = (mp->nspec_outer_inner_core) * NGLL3;
- color_offset_nonpadded_att2 = (mp->nspec_outer_inner_core) * NGLL3 * N_SLS;
+ if( mp->attenuation_3D ){
+ color_offset_nonpadded_att2 = (mp->nspec_outer_inner_core) * NGLL3 * N_SLS;
+ }else{
+ color_offset_nonpadded_att2 = (mp->nspec_outer_inner_core) * 1 * N_SLS;
+ }
color_offset_ispec = mp->nspec_outer_inner_core;
}
@@ -1286,7 +1300,7 @@
//}
Kernel_2_inner_core(nb_blocks_to_compute,mp,
- *deltat,
+ *deltat,
*iphase,
mp->d_ibool_inner_core + color_offset_nonpadded,
mp->d_idoubling_inner_core + color_offset_ispec,
@@ -1337,7 +1351,11 @@
// for no-aligned arrays
color_offset_nonpadded += nb_blocks_to_compute * NGLL3;
// for factor_common array
- color_offset_nonpadded_att2 += nb_blocks_to_compute * NGLL3 * N_SLS;
+ if( mp->attenuation_3D ){
+ color_offset_nonpadded_att2 += nb_blocks_to_compute * NGLL3 * N_SLS;
+ }else{
+ color_offset_nonpadded_att2 += nb_blocks_to_compute * 1 * N_SLS;
+ }
// for array(ispec)
color_offset_ispec += nb_blocks_to_compute;
}
@@ -1347,7 +1365,7 @@
// no mesh coloring: uses atomic updates
Kernel_2_inner_core(num_elements,mp,
- *deltat,
+ *deltat,
*iphase,
mp->d_ibool_inner_core,
mp->d_idoubling_inner_core,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_kernels_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_kernels_cuda.cu 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_kernels_cuda.cu 2012-07-20 21:52:01 UTC (rev 20532)
@@ -100,6 +100,8 @@
}
}
+/* ----------------------------------------------------------------------------------------------- */
+
__device__ void compute_strain_product(realw* prod,
realw eps_trace_over_3,
realw* epsdev,
@@ -150,6 +152,8 @@
}
}
+/* ----------------------------------------------------------------------------------------------- */
+
__global__ void compute_kernels_ani_cudakernel(int* ibool,
realw* epsilondev_xx,
realw* epsilondev_yy,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-07-20 21:52:01 UTC (rev 20532)
@@ -461,9 +461,12 @@
// simulation flags
int save_forward;
int absorbing_conditions;
+
int attenuation;
int attenuation_new;
int use_attenuation_mimic;
+ int attenuation_3D;
+
int compute_and_store_strain;
int anisotropic_3D_mantle;
int gravity;
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-07-20 21:52:01 UTC (rev 20532)
@@ -78,6 +78,7 @@
int* ATTENUATION_f,
int* ATTENUATION_NEW_f,
int* USE_ATTENUATION_MIMIC_f,
+ int* ATTENUATION_3D_VAL_f,
int* COMPUTE_AND_STORE_STRAIN_f,
int* ANISOTROPIC_3D_MANTLE_f,
int* ANISOTROPIC_INNER_CORE_f,
@@ -129,9 +130,12 @@
mp->oceans = *OCEANS_f;
mp->gravity = *GRAVITY_f;
mp->rotation = *ROTATION_f;
+
mp->attenuation = *ATTENUATION_f;
mp->attenuation_new = *ATTENUATION_NEW_f;
mp->use_attenuation_mimic = *USE_ATTENUATION_MIMIC_f;
+ mp->attenuation_3D = *ATTENUATION_3D_VAL_f;
+
mp->compute_and_store_strain = *COMPUTE_AND_STORE_STRAIN_f;
mp->anisotropic_3D_mantle = *ANISOTROPIC_3D_MANTLE_f;
mp->anisotropic_inner_core = *ANISOTROPIC_INNER_CORE_f;
@@ -397,10 +401,16 @@
if( ! mp->attenuation ){ exit_on_cuda_error("prepare_fields_attenuat_device attenuation not properly initialized"); }
// crust_mantle
- R_size1 = N_SLS*NGLL3*mp->NSPEC_CRUST_MANTLE;
- R_size2 = NGLL3*mp->NSPEC_CRUST_MANTLE;
- R_size3 = N_SLS*NGLL3*mp->NSPEC_CRUST_MANTLE;
-
+ if( mp->attenuation_3D ){
+ R_size1 = N_SLS*NGLL3*mp->NSPEC_CRUST_MANTLE;
+ R_size2 = NGLL3*mp->NSPEC_CRUST_MANTLE;
+ R_size3 = N_SLS*NGLL3*mp->NSPEC_CRUST_MANTLE;
+ }else{
+ R_size1 = N_SLS*NGLL3*mp->NSPEC_CRUST_MANTLE;
+ R_size2 = 1*mp->NSPEC_CRUST_MANTLE;
+ R_size3 = N_SLS*1*mp->NSPEC_CRUST_MANTLE;
+ }
+
print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_one_minus_sum_beta_crust_mantle,
R_size2*sizeof(realw)),4430);
print_CUDA_error_if_any(cudaMemcpy(mp->d_one_minus_sum_beta_crust_mantle,one_minus_sum_beta_crust_mantle,
@@ -438,10 +448,16 @@
}
// inner_core
- R_size1 = 5*N_SLS*NGLL3*mp->NSPEC_INNER_CORE;
- R_size2 = NGLL3*mp->NSPEC_INNER_CORE;
- R_size3 = N_SLS*NGLL3*mp->NSPEC_INNER_CORE;
-
+ if( mp->attenuation_3D ){
+ R_size1 = N_SLS*NGLL3*mp->NSPEC_INNER_CORE;
+ R_size2 = NGLL3*mp->NSPEC_INNER_CORE;
+ R_size3 = N_SLS*NGLL3*mp->NSPEC_INNER_CORE;
+ }else{
+ R_size1 = N_SLS*NGLL3*mp->NSPEC_INNER_CORE;
+ R_size2 = 1*mp->NSPEC_INNER_CORE;
+ R_size3 = N_SLS*1*mp->NSPEC_INNER_CORE;
+ }
+
print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_one_minus_sum_beta_inner_core,
R_size2*sizeof(realw)),4430);
print_CUDA_error_if_any(cudaMemcpy(mp->d_one_minus_sum_beta_inner_core,one_minus_sum_beta_inner_core,
@@ -456,7 +472,7 @@
// memory variables
print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_R_xx_inner_core,
- R_size1*sizeof(realw)),4401);
+ R_size1*sizeof(realw)),4401);
print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_R_yy_inner_core,
R_size1*sizeof(realw)),4401);
print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_R_xy_inner_core,
@@ -467,7 +483,7 @@
R_size1*sizeof(realw)),4401);
print_CUDA_error_if_any(cudaMemcpy(mp->d_R_xx_inner_core,R_xx_inner_core,
- R_size1*sizeof(realw),cudaMemcpyHostToDevice),4402);
+ R_size1*sizeof(realw),cudaMemcpyHostToDevice),4402);
print_CUDA_error_if_any(cudaMemcpy(mp->d_R_yy_inner_core,R_yy_inner_core,
R_size1*sizeof(realw),cudaMemcpyHostToDevice),4402);
print_CUDA_error_if_any(cudaMemcpy(mp->d_R_xy_inner_core,R_xy_inner_core,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -269,7 +269,7 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
nspec_ani,nspec_stacey,nspec_att,Qmu_store,tau_e_store,tau_s,T_c_source, &
- size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4), &
+ size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4), &
rho_vp,rho_vs,ACTUALLY_STORE_ARRAYS, &
xigll,yigll,zigll,ispec_is_tiso)
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -262,7 +262,7 @@
ipass)
-! initialize number of layers
+ ! initialize number of layers
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...setting up layers '
@@ -271,7 +271,7 @@
xstore,ystore,zstore,ibool,idoubling, &
NEX_PER_PROC_ETA,is_on_a_slice_edge)
-! creates mesh elements
+ ! creates mesh elements
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...creating mesh elements '
@@ -576,21 +576,22 @@
T_c_source = AM_V%QT_c_source
tau_s(:) = AM_V%Qtau_s(:)
nspec_att = nspec
+
+ if (ATTENUATION_3D) then
+ ! attenuation arrays are fully 3D
+ allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec_att), &
+ tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec_att),stat=ier)
+ else if (.not. ATTENUATION_3D) then
+ ! save some memory in case of 1D attenuation models
+ allocate(Qmu_store(1,1,1,nspec_att), &
+ tau_e_store(N_SLS,1,1,1,nspec_att),stat=ier)
+ endif
else
+ ! allocates dummy size arrays
nspec_att = 1
+ allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec_att), &
+ tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec_att),stat=ier)
end if
-
- if (ATTENUATION) then
- if (ATTENUATION_3D) then
- ! attenuation arrays are fully 3D
- allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec_att), &
- tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec_att),stat=ier)
- else if (.not. ATTENUATION_3D) then
- ! save some memory in case of 1D attenuation models
- allocate(Qmu_store(1,1,1,nspec_att), &
- tau_e_store(N_SLS,1,1,1,nspec_att),stat=ier)
- endif
- endif
if(ier /= 0) stop 'error in allocate 1'
! array with model density
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -42,54 +42,49 @@
implicit none
- integer myrank,iregion_code,ispec,nspec,idoubling
+ integer :: myrank,iregion_code,ispec,nspec,idoubling
- real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) kappahstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) muhstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) eta_anisostore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) dvpstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,nspec) :: kappavstore,kappahstore, &
+ muvstore,muhstore,eta_anisostore,rhostore,dvpstore
- integer nspec_ani
+ integer :: nspec_ani
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_ani) :: &
c11store,c12store,c13store,c14store,c15store,c16store, &
c22store,c23store,c24store,c25store,c26store, &
c33store,c34store,c35store,c36store, &
c44store,c45store,c46store,c55store,c56store,c66store
- integer nspec_stacey
- real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,nspec_stacey),rho_vs(NGLLX,NGLLY,NGLLZ,nspec_stacey)
+ integer :: nspec_stacey
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,nspec_stacey):: rho_vp,rho_vs
double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
- double precision rmin,rmax,RCMB,RICB,R670,RMOHO, &
+ double precision :: rmin,rmax,RCMB,RICB,R670,RMOHO, &
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN
! attenuation values
- integer vx,vy,vz,vnspec
+ integer :: vx,vy,vz,vnspec
double precision, dimension(N_SLS) :: tau_s
double precision, dimension(vx, vy, vz, vnspec) :: Qmu_store
double precision, dimension(N_SLS, vx, vy, vz, vnspec) :: tau_e_store
- double precision T_c_source
+ double precision :: T_c_source
- logical ABSORBING_CONDITIONS
- logical elem_in_crust,elem_in_mantle
+ logical :: ABSORBING_CONDITIONS
+ logical :: elem_in_crust,elem_in_mantle
! local parameters
- double precision xmesh,ymesh,zmesh
+ double precision :: xmesh,ymesh,zmesh
! the 21 coefficients for an anisotropic medium in reduced notation
- double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
+ double precision :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
c34,c35,c36,c44,c45,c46,c55,c56,c66
double precision, dimension(N_SLS) :: tau_e
! local parameters
- double precision rho,dvp
- double precision vpv,vph,vsv,vsh,eta_aniso
- double precision Qkappa,Qmu
- double precision r,r_prem,moho
- integer i,j,k
+ double precision :: rho,dvp
+ double precision :: vpv,vph,vsv,vsh,eta_aniso
+ double precision :: Qkappa,Qmu
+ double precision :: r,r_prem,moho
+ integer :: i,j,k
! loops over all gll points for this spectral element
do k=1,NGLLZ
@@ -312,19 +307,27 @@
endif !CUSTOM_REAL
- if(ATTENUATION .and. ATTENUATION_3D) then
- tau_e_store(:,i,j,k,ispec) = tau_e(:)
- Qmu_store(i,j,k,ispec) = Qmu
+ if( ATTENUATION ) then
+ if( ATTENUATION_3D ) then
+ tau_e_store(:,i,j,k,ispec) = tau_e(:)
+ Qmu_store(i,j,k,ispec) = Qmu
+ else
+ ! store values from mid-point for whole element
+ if( i == NGLLX/2 .and. j == NGLLY/2 .and. k == NGLLZ/2 ) then
+ tau_e_store(:,1,1,1,ispec) = tau_e(:)
+ Qmu_store(1,1,1,ispec) = Qmu
+ endif
+ endif
endif
enddo
enddo
enddo
- if (ATTENUATION .and. .not. ATTENUATION_3D) then
- tau_e_store(:,1,1,1,ispec) = tau_e(:)
- Qmu_store(1,1,1,ispec) = Qmu
- endif
+ !if (ATTENUATION .and. .not. ATTENUATION_3D) then
+ ! tau_e_store(:,1,1,1,ispec) = tau_e(:)
+ ! Qmu_store(1,1,1,ispec) = Qmu
+ !endif
end subroutine get_model
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -356,6 +356,7 @@
!
subroutine model_attenuation_getstored_tau(Qmu_in, T_c_source, tau_s, tau_e, AM_V, AM_S, AS_V)
+
! includes min_period, max_period, and N_SLS
implicit none
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -241,7 +241,7 @@
static_memory_size = static_memory_size + &
6.d0*nglob(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
else
- ! one only keeps one mass matrix for the calculations: rmassz
+ ! one only keeps one mass matrix for the calculations: rmassz
static_memory_size = static_memory_size + &
4.d0*nglob(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -80,7 +80,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
kappavstore,muvstore
- ! variable sized array variables
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
! attenuation
@@ -234,10 +234,12 @@
endif
! precompute terms for attenuation if needed
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
endif
!
@@ -473,7 +475,7 @@
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
- ! variable sized array variables
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
@@ -634,10 +636,12 @@
endif
! precompute terms for attenuation if needed
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
endif
@@ -1214,12 +1218,14 @@
endif
! precompute terms for attenuation if needed
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
- minus_sum_beta = one_minus_sum_beta_use - 1.0
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
- minus_sum_beta = one_minus_sum_beta_use - 1.0
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ minus_sum_beta = one_minus_sum_beta_use - 1.0_CUSTOM_REAL
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ minus_sum_beta = one_minus_sum_beta_use - 1.0_CUSTOM_REAL
+ endif
endif
!
@@ -1572,10 +1578,9 @@
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
-! real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_xx,R_yy,R_xy,R_xz,R_yz
- ! variable sized array variables
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
@@ -1584,7 +1589,6 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: c44store
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: muvstore
-! real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
@@ -1614,84 +1618,84 @@
gammal = gammaval(i_SLS)
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if (ATTENUATION_3D_VAL) then
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
- endif
- else if (.not. ATTENUATION_3D_VAL) then
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
- endif
+ if( ATTENUATION_3D_VAL ) then
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
+ endif
+ else
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
+ endif
endif
! this helps to vectorize the inner most loop
do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ! R_memory(:,i_SLS,i,j,k,ispec) = alphal * R_memory(:,i_SLS,i,j,k,ispec) &
- ! + factor_common_c44_muv(i,j,k) &
- ! *( betal * epsilondev(:,i,j,k,ispec) + gammal * epsilondev_loc(:,i,j,k))
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! R_memory(:,i_SLS,i,j,k,ispec) = alphal * R_memory(:,i_SLS,i,j,k,ispec) &
+ ! + factor_common_c44_muv(i,j,k) &
+ ! *( betal * epsilondev(:,i,j,k,ispec) + gammal * epsilondev_loc(:,i,j,k))
- R_xx(i_SLS,i,j,k,ispec) = alphal * R_xx(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
+ R_xx(i_SLS,i,j,k,ispec) = alphal * R_xx(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
(betal * epsilondev_xx(i,j,k,ispec) + gammal * epsilondev_loc(1,i,j,k))
- R_yy(i_SLS,i,j,k,ispec) = alphal * R_yy(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
+ R_yy(i_SLS,i,j,k,ispec) = alphal * R_yy(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
(betal * epsilondev_yy(i,j,k,ispec) + gammal * epsilondev_loc(2,i,j,k))
- R_xy(i_SLS,i,j,k,ispec) = alphal * R_xy(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
+ R_xy(i_SLS,i,j,k,ispec) = alphal * R_xy(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
(betal * epsilondev_xy(i,j,k,ispec) + gammal * epsilondev_loc(3,i,j,k))
- R_xz(i_SLS,i,j,k,ispec) = alphal * R_xz(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
+ R_xz(i_SLS,i,j,k,ispec) = alphal * R_xz(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
(betal * epsilondev_xz(i,j,k,ispec) + gammal * epsilondev_loc(4,i,j,k))
- R_yz(i_SLS,i,j,k,ispec) = alphal * R_yz(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
+ R_yz(i_SLS,i,j,k,ispec) = alphal * R_yz(i_SLS,i,j,k,ispec) + factor_common_c44_muv(i,j,k) * &
(betal * epsilondev_yz(i,j,k,ispec) + gammal * epsilondev_loc(5,i,j,k))
- enddo
- enddo
+ enddo
+ enddo
enddo
enddo ! i_SLS
#else
! way 1:
do i_SLS = 1,N_SLS
- ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if (ATTENUATION_3D_VAL) then
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
- endif
- else if (.not. ATTENUATION_3D_VAL) then
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
- endif
- endif
+ ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
+ if( ATTENUATION_3D_VAL ) then
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
+ endif
+ else
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
+ endif
+ endif
- ! do i_memory = 1,5
- ! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
- ! + factor_common_c44_muv(:,:,:) &
- ! * (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
- ! enddo
+ ! do i_memory = 1,5
+ ! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ ! + factor_common_c44_muv(:,:,:) &
+ ! * (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ ! enddo
- R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
- R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
- R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
- R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
- R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
enddo ! i_SLS
@@ -1739,12 +1743,10 @@
! element id
integer :: ispec
-! real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
-
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: &
R_xx,R_yy,R_xy,R_xz,R_yz
- ! variable sized array variables
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
@@ -1783,33 +1785,33 @@
gammal = gammaval(i_SLS)
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if (ATTENUATION_3D_VAL) then
- factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
- else if (.not. ATTENUATION_3D_VAL) then
- factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
+ if( ATTENUATION_3D_VAL ) then
+ factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
+ else
+ factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
endif
! this helps to vectorize the inner most loop
do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ! R_memory(:,i_SLS,i,j,k,ispec) = alphal * R_memory(:,i_SLS,i,j,k,ispec) &
- ! + factor_common_use(i,j,k) &
- ! *( betal * epsilondev(:,i,j,k,ispec) + gammal * epsilondev_loc(:,i,j,k))
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! R_memory(:,i_SLS,i,j,k,ispec) = alphal * R_memory(:,i_SLS,i,j,k,ispec) &
+ ! + factor_common_use(i,j,k) &
+ ! *( betal * epsilondev(:,i,j,k,ispec) + gammal * epsilondev_loc(:,i,j,k))
- R_xx(i_SLS,i,j,k,ispec) = alphal * R_xx(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
+ R_xx(i_SLS,i,j,k,ispec) = alphal * R_xx(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
(betal * epsilondev_xx(i,j,k,ispec) + gammal * epsilondev_loc(1,i,j,k))
- R_yy(i_SLS,i,j,k,ispec) = alphal * R_yy(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
+ R_yy(i_SLS,i,j,k,ispec) = alphal * R_yy(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
(betal * epsilondev_yy(i,j,k,ispec) + gammal * epsilondev_loc(2,i,j,k))
- R_xy(i_SLS,i,j,k,ispec) = alphal * R_xy(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
+ R_xy(i_SLS,i,j,k,ispec) = alphal * R_xy(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
(betal * epsilondev_xy(i,j,k,ispec) + gammal * epsilondev_loc(3,i,j,k))
- R_xz(i_SLS,i,j,k,ispec) = alphal * R_xz(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
+ R_xz(i_SLS,i,j,k,ispec) = alphal * R_xz(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
(betal * epsilondev_xz(i,j,k,ispec) + gammal * epsilondev_loc(4,i,j,k))
- R_yz(i_SLS,i,j,k,ispec) = alphal * R_yz(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
+ R_yz(i_SLS,i,j,k,ispec) = alphal * R_yz(i_SLS,i,j,k,ispec) + factor_common_use(i,j,k) * &
(betal * epsilondev_yz(i,j,k,ispec) + gammal * epsilondev_loc(5,i,j,k))
- enddo
- enddo
+ enddo
+ enddo
enddo
enddo ! i_SLS
@@ -1818,10 +1820,10 @@
do i_SLS = 1,N_SLS
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if (ATTENUATION_3D_VAL) then
- factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
- else if (.not. ATTENUATION_3D_VAL) then
- factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
+ if( ATTENUATION_3D_VAL ) then
+ factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
+ else
+ factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
endif
! do i_memory = 1,5
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -79,19 +79,16 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: accel_crust_mantle
- ! variable sized array variables
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
-! real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATT) :: R_xx,R_yy,R_xy,R_xz,R_yz
-! real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
-
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: epsilon_trace_over_3
! [alpha,beta,gamma]val reduced to N_SLS and factor_common to N_SLS*NUM_NODES
@@ -362,12 +359,14 @@
endif
! precompute terms for attenuation if needed
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
endif
- minus_sum_beta = one_minus_sum_beta_use - 1.0
+ minus_sum_beta = one_minus_sum_beta_use - 1.0_CUSTOM_REAL
!
! compute either isotropic or anisotropic elements
@@ -892,19 +891,19 @@
! IMPROVE we should probably use an average value instead
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if (ATTENUATION_3D_VAL) then
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
- endif
- else if (.not. ATTENUATION_3D_VAL) then
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
- endif
- endif
+ if(ATTENUATION_3D_VAL) then
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
+ endif
+ else
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
+ endif
+ endif
! do i_memory = 1,5
! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * &
@@ -914,19 +913,19 @@
! gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
! enddo
- R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
- R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
- R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
- R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
- R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
+ R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
(betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -89,19 +89,16 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: accel_crust_mantle
- ! variable sized array variables
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
-! real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATT) :: R_xx,R_yy,R_xy,R_xz,R_yz
-! real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
-
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: epsilon_trace_over_3
! [alpha,beta,gamma]val reduced to N_SLS and factor_common to N_SLS*NUM_NODES
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -77,19 +77,16 @@
! to allow for optimization of cache access by compiler
! variable lengths for factor_common and one_minus_sum_beta
- ! variable sized array variables
+ ! variable sized array variables
integer vx, vy, vz, vnspec
real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
-! real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATT) :: R_xx,R_yy,R_xy,R_xz,R_yz
-! real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
-
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: epsilon_trace_over_3
! inner/outer element run flag
@@ -346,10 +343,12 @@
endif
! precompute terms for attenuation if needed
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0_CUSTOM_REAL
+ else
+ minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0_CUSTOM_REAL
+ endif
endif
if(ANISOTROPIC_INNER_CORE_VAL) then
@@ -401,10 +400,12 @@
mul = muvstore(i,j,k,ispec)
! use unrelaxed parameters if attenuation
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- mul = mul * one_minus_sum_beta(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- mul = mul * one_minus_sum_beta(1,1,1,ispec)
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ mul = mul * one_minus_sum_beta(i,j,k,ispec)
+ else
+ mul = mul * one_minus_sum_beta(1,1,1,ispec)
+ endif
endif
lambdalplus2mul = kappal + FOUR_THIRDS * mul
@@ -652,12 +653,12 @@
do i_SLS = 1,N_SLS
- ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if (ATTENUATION_3D_VAL) then
- factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
- else if (.not. ATTENUATION_3D_VAL) then
- factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
- endif
+ ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
+ if (ATTENUATION_3D_VAL) then
+ factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
+ else
+ factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
+ endif
! do i_memory = 1,5
! R_memory(i_memory,i_SLS,:,:,:,ispec) = &
@@ -667,20 +668,20 @@
! (betaval(i_SLS) * &
! epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
! enddo
-
- R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
+
+ R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
- R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
+ R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
- R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
+ R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
- R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
+ R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
- R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
+ R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -90,14 +90,11 @@
real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
-! real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATT) :: &
R_xx,R_yy,R_xy,R_xz,R_yz
-! real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
-
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: epsilon_trace_over_3
! inner/outer element run flag
@@ -616,10 +613,12 @@
endif
endif
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0_CUSTOM_REAL
+ else
+ minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0_CUSTOM_REAL
+ endif
endif
if(ANISOTROPIC_INNER_CORE_VAL) then
@@ -669,10 +668,12 @@
mul = muvstore(i,j,k,ispec)
! use unrelaxed parameters if attenuation
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- mul = mul * one_minus_sum_beta(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- mul = mul * one_minus_sum_beta(1,1,1,ispec)
+ if( ATTENUATION_VAL ) then
+ if( ATTENUATION_3D_VAL ) then
+ mul = mul * one_minus_sum_beta(i,j,k,ispec)
+ else
+ mul = mul * one_minus_sum_beta(1,1,1,ispec)
+ endif
endif
lambdalplus2mul = kappal + FOUR_THIRDS * mul
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -27,7 +27,8 @@
subroutine get_attenuation_model_3D_or_1D(myrank, prname, one_minus_sum_beta, &
- factor_common, scale_factor, tau_s, vx, vy, vz, vnspec)
+ factor_common, scale_factor, tau_s, &
+ vx, vy, vz, vnspec)
use specfem_par,only: ATTENUATION_VAL, ATTENUATION_3D_VAL
@@ -35,21 +36,26 @@
include 'constants.h'
- integer myrank,vx,vy,vz,vnspec
- character(len=150) prname
+ integer :: myrank,vx,vy,vz,vnspec
+ character(len=150) :: prname
double precision, dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta, scale_factor
double precision, dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
double precision, dimension(N_SLS) :: tau_s
- integer i,j,k,ispec
-
+ ! local parameters
+ integer :: i,j,k,ispec,ier
double precision, dimension(N_SLS) :: tau_e, fc
- double precision omsb, Q_mu, sf, T_c_source, scale_t
+ double precision :: omsb, Q_mu, sf, T_c_source, scale_t
+
+ ! checks if attenuation is on and anything to do
+ if( .not. ATTENUATION_VAL) return
! All of the following reads use the output parameters as their temporary arrays
! use the filename to determine the actual contents of the read
open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', &
- status='old',action='read',form='unformatted')
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file attenuation.bin')
+
read(27) tau_s
read(27) factor_common
read(27) scale_factor
@@ -63,42 +69,42 @@
T_c_source = 1000.0d0 / T_c_source
T_c_source = T_c_source / scale_t
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- do ispec = 1, vnspec
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- tau_e(:) = factor_common(:,i,j,k,ispec)
- Q_mu = scale_factor(i,j,k,ispec)
+ if( ATTENUATION_3D_VAL ) then
+ do ispec = 1, vnspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ tau_e(:) = factor_common(:,i,j,k,ispec)
+ Q_mu = scale_factor(i,j,k,ispec)
- ! Determine the factor_common and one_minus_sum_beta from tau_s and tau_e
- call get_attenuation_property_values(tau_s, tau_e, fc, omsb)
+ ! Determine the factor_common and one_minus_sum_beta from tau_s and tau_e
+ call get_attenuation_property_values(tau_s, tau_e, fc, omsb)
- factor_common(:,i,j,k,ispec) = fc(:)
- one_minus_sum_beta(i,j,k,ispec) = omsb
+ factor_common(:,i,j,k,ispec) = fc(:)
+ one_minus_sum_beta(i,j,k,ispec) = omsb
- ! Determine the "scale_factor" from tau_s, tau_e, central source frequency, and Q
- call get_attenuation_scale_factor(myrank, T_c_source, tau_e, tau_s, Q_mu, sf)
- scale_factor(i,j,k,ispec) = sf
- enddo
- enddo
+ ! Determine the "scale_factor" from tau_s, tau_e, central source frequency, and Q
+ call get_attenuation_scale_factor(myrank, T_c_source, tau_e, tau_s, Q_mu, sf)
+ scale_factor(i,j,k,ispec) = sf
+ enddo
enddo
- enddo
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- do ispec = 1, vnspec
- tau_e(:) = factor_common(:,1,1,1,ispec)
- Q_mu = scale_factor(1,1,1,ispec)
-
- ! Determine the factor_common and one_minus_sum_beta from tau_s and tau_e
- call get_attenuation_property_values(tau_s, tau_e, fc, omsb)
-
- factor_common(:,1,1,1,ispec) = fc(:)
- one_minus_sum_beta(1,1,1,ispec) = omsb
-
- ! Determine the "scale_factor" from tau_s, tau_e, central source frequency, and Q
- call get_attenuation_scale_factor(myrank, T_c_source, tau_e, tau_s, Q_mu, sf)
- scale_factor(1,1,1,ispec) = sf
- enddo
+ enddo
+ enddo
+ else
+ do ispec = 1, vnspec
+ tau_e(:) = factor_common(:,1,1,1,ispec)
+ Q_mu = scale_factor(1,1,1,ispec)
+
+ ! Determine the factor_common and one_minus_sum_beta from tau_s and tau_e
+ call get_attenuation_property_values(tau_s, tau_e, fc, omsb)
+
+ factor_common(:,1,1,1,ispec) = fc(:)
+ one_minus_sum_beta(1,1,1,ispec) = omsb
+
+ ! Determine the "scale_factor" from tau_s, tau_e, central source frequency, and Q
+ call get_attenuation_scale_factor(myrank, T_c_source, tau_e, tau_s, Q_mu, sf)
+ scale_factor(1,1,1,ispec) = sf
+ enddo
endif
end subroutine get_attenuation_model_3D_or_1D
@@ -214,14 +220,14 @@
double precision, dimension(N_SLS) :: tauinv
- tauinv(:) = - 1.0 / tau_s(:)
+ tauinv(:) = - 1.d0 / tau_s(:)
- alphaval(:) = 1 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2. + &
- deltat**3*tauinv(:)**3 / 6. + deltat**4*tauinv(:)**4 / 24.
- betaval(:) = deltat / 2. + deltat**2*tauinv(:) / 3. &
- + deltat**3*tauinv(:)**2 / 8. + deltat**4*tauinv(:)**3 / 24.
- gammaval(:) = deltat / 2. + deltat**2*tauinv(:) / 6. &
- + deltat**3*tauinv(:)**2 / 24.0
+ alphaval(:) = 1.d0 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2.d0 + &
+ deltat**3*tauinv(:)**3 / 6.d0 + deltat**4*tauinv(:)**4 / 24.d0
+ betaval(:) = deltat / 2.d0 + deltat**2*tauinv(:) / 3.d0 &
+ + deltat**3*tauinv(:)**2 / 8.d0 + deltat**4*tauinv(:)**3 / 24.d0
+ gammaval(:) = deltat / 2.d0 + deltat**2*tauinv(:) / 6.d0 &
+ + deltat**3*tauinv(:)**2 / 24.d0
end subroutine get_attenuation_memory_values
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -64,7 +64,7 @@
call prepare_timerun_gravity()
! precomputes attenuation factors
- if(ATTENUATION_VAL) call prepare_timerun_attenuation()
+ call prepare_timerun_attenuation()
! initializes arrays
call prepare_timerun_init_wavefield()
@@ -682,6 +682,9 @@
integer :: ispec,i,j,k,ier
character(len=150) :: prnamel
+ ! checks if attenuation is on and anything to do
+ if( .not. ATTENUATION_VAL ) return
+
! get and store PREM attenuation model
if(myrank == 0 ) then
write(IMAIN,*) "preparing attenuation."
@@ -705,15 +708,29 @@
! CRUST_MANTLE ATTENUATION
call create_name_database(prnamel, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
+
+ ! initializes
+ omsb_crust_mantle_dble = 0.d0
+ factor_common_crust_mantle_dble = 0.d0
+ factor_scale_crust_mantle_dble = 0.d0
+ tau_sigma_dble = 0.d0
+
call get_attenuation_model_3D_or_1D(myrank, prnamel, omsb_crust_mantle_dble, &
factor_common_crust_mantle_dble,factor_scale_crust_mantle_dble,tau_sigma_dble, &
- ATT1,ATT2,ATT3,ATT4,NSPEC_CRUST_MANTLE)
+ ATT1,ATT2,ATT3,ATT4)
! INNER_CORE ATTENUATION
call create_name_database(prnamel, myrank, IREGION_INNER_CORE, LOCAL_PATH)
+
+ ! initializes
+ omsb_inner_core_dble = 0.d0
+ factor_common_inner_core_dble = 0.d0
+ factor_scale_inner_core_dble = 0.d0
+ tau_sigma_dble = 0.d0
+
call get_attenuation_model_3D_or_1D(myrank, prnamel, omsb_inner_core_dble, &
factor_common_inner_core_dble,factor_scale_inner_core_dble,tau_sigma_dble, &
- ATT1,ATT2,ATT3,ATT5,NSPEC_INNER_CORE)
+ ATT1,ATT2,ATT3,ATT5)
if(CUSTOM_REAL == SIZE_REAL) then
factor_scale_crust_mantle = sngl(factor_scale_crust_mantle_dble)
@@ -748,67 +765,64 @@
! W. H. Freeman, (1980), second edition, sections 5.5 and 5.5.2, eq. (5.81) p. 170
! rescale in crust and mantle
-
do ispec = 1,NSPEC_CRUST_MANTLE
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- scale_factor = factor_scale_crust_mantle(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- scale_factor = factor_scale_crust_mantle(1,1,1,ispec)
- endif
-
- if(ANISOTROPIC_3D_MANTLE_VAL) then
- scale_factor_minus_one = scale_factor - 1.
- mul = c44store_crust_mantle(i,j,k,ispec)
- c11store_crust_mantle(i,j,k,ispec) = c11store_crust_mantle(i,j,k,ispec) &
- + FOUR_THIRDS * scale_factor_minus_one * mul
- c12store_crust_mantle(i,j,k,ispec) = c12store_crust_mantle(i,j,k,ispec) &
- - TWO_THIRDS * scale_factor_minus_one * mul
- c13store_crust_mantle(i,j,k,ispec) = c13store_crust_mantle(i,j,k,ispec) &
- - TWO_THIRDS * scale_factor_minus_one * mul
- c22store_crust_mantle(i,j,k,ispec) = c22store_crust_mantle(i,j,k,ispec) &
- + FOUR_THIRDS * scale_factor_minus_one * mul
- c23store_crust_mantle(i,j,k,ispec) = c23store_crust_mantle(i,j,k,ispec) &
- - TWO_THIRDS * scale_factor_minus_one * mul
- c33store_crust_mantle(i,j,k,ispec) = c33store_crust_mantle(i,j,k,ispec) &
- + FOUR_THIRDS * scale_factor_minus_one * mul
- c44store_crust_mantle(i,j,k,ispec) = c44store_crust_mantle(i,j,k,ispec) &
- + scale_factor_minus_one * mul
- c55store_crust_mantle(i,j,k,ispec) = c55store_crust_mantle(i,j,k,ispec) &
- + scale_factor_minus_one * mul
- c66store_crust_mantle(i,j,k,ispec) = c66store_crust_mantle(i,j,k,ispec) &
- + scale_factor_minus_one * mul
- else
- if(MOVIE_VOLUME .and. SIMULATION_TYPE==3) then
- ! store the original value of \mu to comput \mu*\eps
- muvstore_crust_mantle_3dmovie(i,j,k,ispec)=muvstore_crust_mantle(i,j,k,ispec)
- endif
- muvstore_crust_mantle(i,j,k,ispec) = muvstore_crust_mantle(i,j,k,ispec) * scale_factor
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if( ATTENUATION_3D_VAL ) then
+ scale_factor = factor_scale_crust_mantle(i,j,k,ispec)
+ else
+ scale_factor = factor_scale_crust_mantle(1,1,1,ispec)
+ endif
- ! scales transverse isotropic values for mu_h
- if( ispec_is_tiso_crust_mantle(ispec) ) then
- muhstore_crust_mantle(i,j,k,ispec) = muhstore_crust_mantle(i,j,k,ispec) * scale_factor
- endif
- endif
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ scale_factor_minus_one = scale_factor - 1.d0
+ mul = c44store_crust_mantle(i,j,k,ispec)
+ c11store_crust_mantle(i,j,k,ispec) = c11store_crust_mantle(i,j,k,ispec) &
+ + FOUR_THIRDS * scale_factor_minus_one * mul
+ c12store_crust_mantle(i,j,k,ispec) = c12store_crust_mantle(i,j,k,ispec) &
+ - TWO_THIRDS * scale_factor_minus_one * mul
+ c13store_crust_mantle(i,j,k,ispec) = c13store_crust_mantle(i,j,k,ispec) &
+ - TWO_THIRDS * scale_factor_minus_one * mul
+ c22store_crust_mantle(i,j,k,ispec) = c22store_crust_mantle(i,j,k,ispec) &
+ + FOUR_THIRDS * scale_factor_minus_one * mul
+ c23store_crust_mantle(i,j,k,ispec) = c23store_crust_mantle(i,j,k,ispec) &
+ - TWO_THIRDS * scale_factor_minus_one * mul
+ c33store_crust_mantle(i,j,k,ispec) = c33store_crust_mantle(i,j,k,ispec) &
+ + FOUR_THIRDS * scale_factor_minus_one * mul
+ c44store_crust_mantle(i,j,k,ispec) = c44store_crust_mantle(i,j,k,ispec) &
+ + scale_factor_minus_one * mul
+ c55store_crust_mantle(i,j,k,ispec) = c55store_crust_mantle(i,j,k,ispec) &
+ + scale_factor_minus_one * mul
+ c66store_crust_mantle(i,j,k,ispec) = c66store_crust_mantle(i,j,k,ispec) &
+ + scale_factor_minus_one * mul
+ else
+ if(MOVIE_VOLUME .and. SIMULATION_TYPE==3) then
+ ! store the original value of \mu to comput \mu*\eps
+ muvstore_crust_mantle_3dmovie(i,j,k,ispec)=muvstore_crust_mantle(i,j,k,ispec)
+ endif
+ muvstore_crust_mantle(i,j,k,ispec) = muvstore_crust_mantle(i,j,k,ispec) * scale_factor
- enddo
+ ! scales transverse isotropic values for mu_h
+ if( ispec_is_tiso_crust_mantle(ispec) ) then
+ muhstore_crust_mantle(i,j,k,ispec) = muhstore_crust_mantle(i,j,k,ispec) * scale_factor
+ endif
+ endif
enddo
- enddo
+ enddo
+ enddo
enddo ! END DO CRUST MANTLE
! rescale in inner core
-
do ispec = 1,NSPEC_INNER_CORE
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- scale_factor_minus_one = factor_scale_inner_core(i,j,k,ispec) - 1.0
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- scale_factor_minus_one = factor_scale_inner_core(1,1,1,ispec) - 1.0
- endif
+ if( ATTENUATION_3D_VAL ) then
+ scale_factor_minus_one = factor_scale_inner_core(i,j,k,ispec) - 1.d0
+ else
+ scale_factor_minus_one = factor_scale_inner_core(1,1,1,ispec) - 1.d0
+ endif
if(ANISOTROPIC_INNER_CORE_VAL) then
mul = muvstore_inner_core(i,j,k,ispec)
@@ -824,12 +838,11 @@
+ scale_factor_minus_one * mul
endif
- if (ATTENUATION_VAL .and. ATTENUATION_3D_VAL) then
- muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(i,j,k,ispec)
- else if (ATTENUATION_VAL .and. .not. ATTENUATION_3D_VAL) then
- muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(1,1,1,ispec)
+ if( ATTENUATION_3D_VAL ) then
+ muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(i,j,k,ispec)
+ else
+ muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(1,1,1,ispec)
endif
-
enddo
enddo
enddo
@@ -860,6 +873,7 @@
endif
endif
+ ! synchronizes processes
call sync_all()
end subroutine prepare_timerun_attenuation
@@ -1134,7 +1148,8 @@
SIMULATION_TYPE,NOISE_TOMOGRAPHY, &
SAVE_FORWARD,ABSORBING_CONDITIONS, &
OCEANS_VAL,GRAVITY_VAL,ROTATION_VAL, &
- ATTENUATION_VAL,ATTENUATION_NEW_VAL,USE_ATTENUATION_MIMIC, &
+ ATTENUATION_VAL,ATTENUATION_NEW_VAL, &
+ USE_ATTENUATION_MIMIC,ATTENUATION_3D_VAL, &
COMPUTE_AND_STORE_STRAIN, &
ANISOTROPIC_3D_MANTLE_VAL,ANISOTROPIC_INNER_CORE_VAL, &
SAVE_BOUNDARY_MESH, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2012-07-20 15:52:51 UTC (rev 20531)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2012-07-20 21:52:01 UTC (rev 20532)
@@ -250,13 +250,13 @@
read(55) b_R_xy_crust_mantle
read(55) b_R_xz_crust_mantle
read(55) b_R_yz_crust_mantle
-
+
read(55) b_R_xx_inner_core
read(55) b_R_yy_inner_core
read(55) b_R_xy_inner_core
read(55) b_R_xz_inner_core
read(55) b_R_yz_inner_core
-
+
! note: for kernel simulations (SIMULATION_TYPE == 3), attenuation is
! only mimicking effects on phase shifts, but not on amplitudes.
! flag USE_ATTENUATION_MIMIC will have to be set to true in this case.
@@ -265,7 +265,7 @@
! therefore no need to transfer arrays onto GPU
!if(GPU_MODE) then
!endif
-
+
endif
close(55)
More information about the CIG-COMMITS
mailing list