[cig-commits] r20375 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: EXAMPLES/regional_Greece_small/DATA src/cuda src/meshfem3D src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Fri Jun 15 04:07:04 PDT 2012
Author: joseph.charles
Date: 2012-06-15 04:07:03 -0700 (Fri, 15 Jun 2012)
New Revision: 20375
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
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/it_update_displacement_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_mass_matrices.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/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.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_acoustic.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_elastic.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/compute_stacey_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.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_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
Log:
fixes Stacey boundary conditions bug by adding a C*deltat/2 contribution to the mass matrices on Stacey edges; updates routines for attenuation debugging
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file 2012-06-15 11:07:03 UTC (rev 20375)
@@ -45,7 +45,7 @@
ATTENUATION_NEW = .false.
# absorbing boundary conditions for a regional simulation
-ABSORBING_CONDITIONS = .true.
+ABSORBING_CONDITIONS = .false.
# record length in minutes
RECORD_LENGTH_IN_MINUTES = 2.5d0
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-06-15 11:07:03 UTC (rev 20375)
@@ -560,7 +560,9 @@
__global__ void compute_coupling_ocean_cuda_kernel(realw* accel_crust_mantle,
- realw* rmass_crust_mantle,
+ realw* rmassx_crust_mantle,
+ realw* rmassy_crust_mantle,
+ realw* rmassz_crust_mantle,
realw* rmass_ocean_load,
realw* normal_top_crust_mantle,
int* ibool_crust_mantle,
@@ -574,7 +576,8 @@
int k,iglob,ispec;
realw nx,ny,nz;
- realw force_normal_comp,additional_term;
+ realw force_normal_comp;
+ realw additional_term_x,additional_term_y,additional_term_z;
// for surfaces elements exactly at the top of the crust mantle (ocean bottom)
if( iface < NSPEC2D_TOP_CM ){
@@ -606,16 +609,18 @@
// make updated component of right-hand side
// we divide by rmass() which is 1 / M
// we use the total force which includes the Coriolis term above
- force_normal_comp = ( accel_crust_mantle[iglob*3]*nx
- + accel_crust_mantle[iglob*3+1]*ny
- + accel_crust_mantle[iglob*3+2]*nz ) / rmass_crust_mantle[iglob];
+ force_normal_comp = accel_crust_mantle[iglob*3]*nx / rmassx_crust_mantle[iglob]
+ + accel_crust_mantle[iglob*3+1]*ny / rmassy_crust_mantle[iglob]
+ + accel_crust_mantle[iglob*3+2]*nz / rmassz_crust_mantle[iglob];
- additional_term = (rmass_ocean_load[iglob] - rmass_crust_mantle[iglob]) * force_normal_comp;
+ additional_term_x = (rmass_ocean_load[iglob] - rmassx_crust_mantle[iglob]) * force_normal_comp;
+ additional_term_y = (rmass_ocean_load[iglob] - rmassy_crust_mantle[iglob]) * force_normal_comp;
+ additional_term_z = (rmass_ocean_load[iglob] - rmassz_crust_mantle[iglob]) * force_normal_comp;
// probably wouldn't need atomicAdd anymore, but just to be sure...
- atomicAdd(&accel_crust_mantle[iglob*3], + additional_term * nx);
- atomicAdd(&accel_crust_mantle[iglob*3+1], + additional_term * ny);
- atomicAdd(&accel_crust_mantle[iglob*3+2], + additional_term * nz);
+ atomicAdd(&accel_crust_mantle[iglob*3], + additional_term_x * nx);
+ atomicAdd(&accel_crust_mantle[iglob*3+1], + additional_term_y * ny);
+ atomicAdd(&accel_crust_mantle[iglob*3+2], + additional_term_z * nz);
}
}
}
@@ -624,7 +629,8 @@
extern "C"
void FC_FUNC_(compute_coupling_ocean_cuda,
- COMPUTE_COUPLING_OCEAN_CUDA)(long* Mesh_pointer_f) {
+ COMPUTE_COUPLING_OCEAN_CUDA)(long* Mesh_pointer_f,
+ int* NCHUNKS_VAL) {
TRACE("compute_coupling_ocean_cuda");
@@ -648,29 +654,64 @@
exit_on_cuda_error("before kernel compute_coupling_ocean_cuda");
#endif
- compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_accel_crust_mantle,
- mp->d_rmass_crust_mantle,
- mp->d_rmass_ocean_load,
- mp->d_normal_top_crust_mantle,
- mp->d_ibool_crust_mantle,
- mp->d_ibelm_top_crust_mantle,
- mp->d_updated_dof_ocean_load,
- mp->nspec2D_top_crust_mantle);
+ if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
+ compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_accel_crust_mantle,
+ mp->d_rmassx_crust_mantle,
+ mp->d_rmassy_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmass_ocean_load,
+ mp->d_normal_top_crust_mantle,
+ mp->d_ibool_crust_mantle,
+ mp->d_ibelm_top_crust_mantle,
+ mp->d_updated_dof_ocean_load,
+ mp->nspec2D_top_crust_mantle);
- // for backward/reconstructed potentials
- if( mp->simulation_type == 3 ) {
- // re-initializes array
- print_CUDA_error_if_any(cudaMemset(mp->d_updated_dof_ocean_load,0,
- sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),88502);
+ // for backward/reconstructed potentials
+ if( mp->simulation_type == 3 ) {
+ // re-initializes array
+ print_CUDA_error_if_any(cudaMemset(mp->d_updated_dof_ocean_load,0,
+ sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),88502);
- compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
- mp->d_rmass_crust_mantle,
- mp->d_rmass_ocean_load,
- mp->d_normal_top_crust_mantle,
- mp->d_ibool_crust_mantle,
- mp->d_ibelm_top_crust_mantle,
- mp->d_updated_dof_ocean_load,
- mp->nspec2D_top_crust_mantle);
+ compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
+ mp->d_rmassx_crust_mantle,
+ mp->d_rmassy_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmass_ocean_load,
+ mp->d_normal_top_crust_mantle,
+ mp->d_ibool_crust_mantle,
+ mp->d_ibelm_top_crust_mantle,
+ mp->d_updated_dof_ocean_load,
+ mp->nspec2D_top_crust_mantle);
+ }
+ }else{
+ compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_accel_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmass_ocean_load,
+ mp->d_normal_top_crust_mantle,
+ mp->d_ibool_crust_mantle,
+ mp->d_ibelm_top_crust_mantle,
+ mp->d_updated_dof_ocean_load,
+ mp->nspec2D_top_crust_mantle);
+
+ // for backward/reconstructed potentials
+ if( mp->simulation_type == 3 ) {
+ // re-initializes array
+ print_CUDA_error_if_any(cudaMemset(mp->d_updated_dof_ocean_load,0,
+ sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),88502);
+
+ compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmass_ocean_load,
+ mp->d_normal_top_crust_mantle,
+ mp->d_ibool_crust_mantle,
+ mp->d_ibelm_top_crust_mantle,
+ mp->d_updated_dof_ocean_load,
+ mp->nspec2D_top_crust_mantle);
+ }
}
#ifdef ENABLE_VERY_SLOW_ERROR_CHECKING
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu 2012-06-15 11:07:03 UTC (rev 20375)
@@ -1101,21 +1101,21 @@
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++) {
+ 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;
@@ -1134,7 +1134,7 @@
tempy3l_att += s_dummyy_loc_att[offset]*hp3;
tempz3l_att += s_dummyz_loc_att[offset]*hp3;
- }
+ }
}
#else
@@ -1193,61 +1193,61 @@
+ 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
@@ -1286,60 +1286,60 @@
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
@@ -1556,9 +1556,9 @@
}
#ifdef USE_TEXTURES
- d_accel[iglob] = tex1Dfetch(tex_accel, iglob) + sum_terms1);
- d_accel[iglob + NGLOB] = tex1Dfetch(tex_accel, iglob + NGLOB) + sum_terms2);
- d_accel[iglob + 2*NGLOB] = tex1Dfetch(tex_accel, iglob + 2*NGLOB) + sum_terms3);
+ d_accel[iglob] = tex1Dfetch(tex_accel, iglob) + sum_terms1 ;
+ d_accel[iglob + NGLOB] = tex1Dfetch(tex_accel, iglob + NGLOB) + sum_terms2 ;
+ d_accel[iglob + 2*NGLOB] = tex1Dfetch(tex_accel, iglob + 2*NGLOB) + sum_terms3 ;
#else
/* OLD/To be implemented version that uses coloring to get around race condition. About 1.6x faster */
@@ -1623,8 +1623,8 @@
d_accel[iglob + NGLOB] -= 0.00000001f;
d_accel[iglob + 2*NGLOB] -= 0.00000001f;
#endif // of #ifndef MAKE_KERNEL2_BECOME_STUPID_FOR_TESTS
+
}
-
/* ----------------------------------------------------------------------------------------------- */
void Kernel_2_crust_mantle(int nb_blocks_to_compute,Mesh* mp,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu 2012-06-15 11:07:03 UTC (rev 20375)
@@ -297,15 +297,17 @@
__global__ void kernel_3_cuda_device(realw* veloc,
realw* accel, int size,
realw deltatover2,
- realw* rmass) {
+ realw* rmassx,
+ realw* rmassy,
+ realw* rmassz) {
int id = threadIdx.x + blockIdx.x*blockDim.x + blockIdx.y*gridDim.x*blockDim.x;
/* because of block and grid sizing problems, there is a small */
/* amount of buffer at the end of the calculation */
if(id < size) {
- accel[3*id] = accel[3*id]*rmass[id];
- accel[3*id+1] = accel[3*id+1]*rmass[id];
- accel[3*id+2] = accel[3*id+2]*rmass[id];
+ accel[3*id] = accel[3*id]*rmassx[id];
+ accel[3*id+1] = accel[3*id+1]*rmassy[id];
+ accel[3*id+2] = accel[3*id+2]*rmassz[id];
veloc[3*id] = veloc[3*id] + deltatover2*accel[3*id];
veloc[3*id+1] = veloc[3*id+1] + deltatover2*accel[3*id+1];
@@ -317,15 +319,17 @@
__global__ void kernel_3_accel_cuda_device(realw* accel,
int size,
- realw* rmass) {
+ realw* rmassx,
+ realw* rmassy,
+ realw* rmassz) {
int id = threadIdx.x + blockIdx.x*blockDim.x + blockIdx.y*gridDim.x*blockDim.x;
/* because of block and grid sizing problems, there is a small */
/* amount of buffer at the end of the calculation */
if(id < size) {
- accel[3*id] = accel[3*id]*rmass[id];
- accel[3*id+1] = accel[3*id+1]*rmass[id];
- accel[3*id+2] = accel[3*id+2]*rmass[id];
+ accel[3*id] = accel[3*id]*rmassx[id];
+ accel[3*id+1] = accel[3*id+1]*rmassy[id];
+ accel[3*id+2] = accel[3*id+2]*rmassz[id];
}
}
@@ -354,7 +358,8 @@
realw* deltatover2_F,
int* SIMULATION_TYPE_f,
realw* b_deltatover2_F,
- int* OCEANS) {
+ int* OCEANS,
+ int* NCHUNKS_VAL) {
TRACE("kernel_3_a_cuda");
Mesh* mp = (Mesh*)(*Mesh_pointer); // get Mesh from fortran integer wrapper
@@ -380,27 +385,76 @@
// check whether we can update accel and veloc, or only accel at this point
if( *OCEANS == 0 ){
// updates both, accel and veloc
- kernel_3_cuda_device<<< grid, threads>>>(mp->d_veloc_crust_mantle,
- mp->d_accel_crust_mantle,
- mp->NGLOB_CRUST_MANTLE,
- deltatover2, mp->d_rmass_crust_mantle);
- if(SIMULATION_TYPE == 3) {
- kernel_3_cuda_device<<< grid, threads>>>(mp->d_b_veloc_crust_mantle,
- mp->d_b_accel_crust_mantle,
- mp->NGLOB_CRUST_MANTLE,
- b_deltatover2,mp->d_rmass_crust_mantle);
+ if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
+ kernel_3_cuda_device<<< grid, threads>>>(mp->d_veloc_crust_mantle,
+ mp->d_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ deltatover2,
+ mp->d_rmassx_crust_mantle,
+ mp->d_rmassy_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+
+ if(SIMULATION_TYPE == 3){
+ kernel_3_cuda_device<<< grid, threads>>>(mp->d_b_veloc_crust_mantle,
+ mp->d_b_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ b_deltatover2,
+ mp->d_rmassx_crust_mantle,
+ mp->d_rmassy_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+ }
+ }else{
+ kernel_3_cuda_device<<< grid, threads>>>(mp->d_veloc_crust_mantle,
+ mp->d_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ deltatover2,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+
+ if(SIMULATION_TYPE == 3){
+ kernel_3_cuda_device<<< grid, threads>>>(mp->d_b_veloc_crust_mantle,
+ mp->d_b_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ b_deltatover2,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+ }
}
+
}else{
// updates only accel
- kernel_3_accel_cuda_device<<< grid, threads>>>(mp->d_accel_crust_mantle,
- mp->NGLOB_CRUST_MANTLE,
- mp->d_rmass_crust_mantle);
- if(SIMULATION_TYPE == 3) {
- kernel_3_accel_cuda_device<<< grid, threads>>>(mp->d_b_accel_crust_mantle,
- mp->NGLOB_CRUST_MANTLE,
- mp->d_rmass_crust_mantle);
+ if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
+ kernel_3_accel_cuda_device<<< grid, threads>>>(mp->d_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ mp->d_rmassx_crust_mantle,
+ mp->d_rmassy_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+
+ if(SIMULATION_TYPE == 3) {
+ kernel_3_accel_cuda_device<<< grid, threads>>>(mp->d_b_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ mp->d_rmassx_crust_mantle,
+ mp->d_rmassy_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+ }
+ }else{
+ kernel_3_accel_cuda_device<<< grid, threads>>>(mp->d_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+
+ if(SIMULATION_TYPE == 3) {
+ kernel_3_accel_cuda_device<<< grid, threads>>>(mp->d_b_accel_crust_mantle,
+ mp->NGLOB_CRUST_MANTLE,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle,
+ mp->d_rmassz_crust_mantle);
+ }
}
}
@@ -472,13 +526,19 @@
kernel_3_cuda_device<<< grid, threads>>>(mp->d_veloc_inner_core,
mp->d_accel_inner_core,
mp->NGLOB_INNER_CORE,
- deltatover2, mp->d_rmass_inner_core);
+ deltatover2,
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core);
if(SIMULATION_TYPE == 3) {
kernel_3_cuda_device<<< grid, threads>>>(mp->d_b_veloc_inner_core,
mp->d_b_accel_inner_core,
mp->NGLOB_INNER_CORE,
- b_deltatover2,mp->d_rmass_inner_core);
+ b_deltatover2,
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core);
}
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-06-15 11:07:03 UTC (rev 20375)
@@ -181,7 +181,9 @@
realw* d_kappavstore_crust_mantle; realw* d_muvstore_crust_mantle;
realw* d_kappahstore_crust_mantle; realw* d_muhstore_crust_mantle;
realw* d_eta_anisostore_crust_mantle;
- realw* d_rmass_crust_mantle;
+ realw* d_rmassx_crust_mantle;
+ realw* d_rmassy_crust_mantle;
+ realw* d_rmassz_crust_mantle;
// global indexing
int* d_ibool_crust_mantle;
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-06-15 11:07:03 UTC (rev 20375)
@@ -1345,7 +1345,7 @@
extern "C"
void FC_FUNC_(prepare_crust_mantle_device,
- PREPARE_CRUST_MANTLE_DEVICE)(long* Mesh_pointer_f,
+ PREPARE_CRUST_MANTLE_DEVICE)(long* Mesh_pointer_f,
realw* h_xix, realw* h_xiy, realw* h_xiz,
realw* h_etax, realw* h_etay, realw* h_etaz,
realw* h_gammax, realw* h_gammay, realw* h_gammaz,
@@ -1353,7 +1353,9 @@
realw* h_kappav, realw* h_muv,
realw* h_kappah, realw* h_muh,
realw* h_eta_aniso,
- realw* h_rmass,
+ realw* h_rmassx,
+ realw* h_rmassy,
+ realw* h_rmassz,
realw* h_normal_top_crust_mantle,
int* h_ibelm_top_crust_mantle,
int* h_ibelm_bottom_crust_mantle,
@@ -1372,7 +1374,8 @@
int* nspec_outer,
int* nspec_inner,
int* NSPEC2D_TOP_CM,
- int* NSPEC2D_BOTTOM_CM
+ int* NSPEC2D_BOTTOM_CM,
+ int* NCHUNKS_VAL
) {
TRACE("prepare_crust_mantle_device");
@@ -1610,9 +1613,18 @@
print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_b_accel_crust_mantle),sizeof(realw)*size),4003);
}
- // mass matrix
- print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_rmass_crust_mantle),sizeof(realw)*size_glob),2005);
- print_CUDA_error_if_any(cudaMemcpy(mp->d_rmass_crust_mantle,h_rmass,
+ // mass matrices
+ if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_rmassx_crust_mantle),sizeof(realw)*size_glob),2005);
+ print_CUDA_error_if_any(cudaMemcpy(mp->d_rmassx_crust_mantle,h_rmassx,
+ sizeof(realw)*size_glob,cudaMemcpyHostToDevice),2100);
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_rmassy_crust_mantle),sizeof(realw)*size_glob),2005);
+ print_CUDA_error_if_any(cudaMemcpy(mp->d_rmassy_crust_mantle,h_rmassy,
+ sizeof(realw)*size_glob,cudaMemcpyHostToDevice),2100);
+ }
+
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_rmassz_crust_mantle),sizeof(realw)*size_glob),2005);
+ print_CUDA_error_if_any(cudaMemcpy(mp->d_rmassz_crust_mantle,h_rmassz,
sizeof(realw)*size_glob,cudaMemcpyHostToDevice),2100);
// kernels
@@ -1669,20 +1681,20 @@
realw* h_gammax, realw* h_gammay, realw* h_gammaz,
realw* h_rho, realw* h_kappav,
realw* h_rmass,
- realw* h_normal_top_outer_core,
- realw* h_normal_bottom_outer_core,
- realw* h_jacobian2D_top_outer_core,
- realw* h_jacobian2D_bottom_outer_core,
- int* h_ibelm_top_outer_core,
- int* h_ibelm_bottom_outer_core,
+ realw* h_normal_top_outer_core,
+ realw* h_normal_bottom_outer_core,
+ realw* h_jacobian2D_top_outer_core,
+ realw* h_jacobian2D_bottom_outer_core,
+ int* h_ibelm_top_outer_core,
+ int* h_ibelm_bottom_outer_core,
int* h_ibool,
realw* h_xstore, realw* h_ystore, realw* h_zstore,
int* num_phase_ispec,
int* phase_ispec_inner,
int* nspec_outer,
int* nspec_inner,
- int* NSPEC2D_TOP_OC,
- int* NSPEC2D_BOTTOM_OC
+ int* NSPEC2D_TOP_OC,
+ int* NSPEC2D_BOTTOM_OC
) {
TRACE("prepare_outer_core_device");
@@ -1836,23 +1848,23 @@
extern "C"
void FC_FUNC_(prepare_inner_core_device,
PREPARE_INNER_CORE_DEVICE)(long* Mesh_pointer_f,
- realw* h_xix, realw* h_xiy, realw* h_xiz,
- realw* h_etax, realw* h_etay, realw* h_etaz,
- realw* h_gammax, realw* h_gammay, realw* h_gammaz,
- realw* h_rho, realw* h_kappav, realw* h_muv,
- realw* h_rmass,
- int* h_ibelm_top_inner_core,
- int* h_ibool,
- realw* h_xstore, realw* h_ystore, realw* h_zstore,
- realw *c11store,realw *c12store,realw *c13store,
- realw *c33store,realw *c44store,
- int* h_idoubling_inner_core,
- int* num_phase_ispec,
- int* phase_ispec_inner,
- int* nspec_outer,
- int* nspec_inner,
- int* NSPEC2D_TOP_IC) {
-
+ realw* h_xix, realw* h_xiy, realw* h_xiz,
+ realw* h_etax, realw* h_etay, realw* h_etaz,
+ realw* h_gammax, realw* h_gammay, realw* h_gammaz,
+ realw* h_rho, realw* h_kappav, realw* h_muv,
+ realw* h_rmass,
+ int* h_ibelm_top_inner_core,
+ int* h_ibool,
+ realw* h_xstore, realw* h_ystore, realw* h_zstore,
+ realw *c11store,realw *c12store,realw *c13store,
+ realw *c33store,realw *c44store,
+ int* h_idoubling_inner_core,
+ int* num_phase_ispec,
+ int* phase_ispec_inner,
+ int* nspec_outer,
+ int* nspec_inner,
+ int* NSPEC2D_TOP_IC) {
+
TRACE("prepare_inner_core_device");
Mesh* mp = (Mesh*)(*Mesh_pointer_f);
@@ -2419,7 +2431,8 @@
extern "C"
void FC_FUNC_(prepare_cleanup_device,
- PREPARE_CLEANUP_DEVICE)(long* Mesh_pointer_f) {
+ PREPARE_CLEANUP_DEVICE)(long* Mesh_pointer_f,
+ int* NCHUNKS_VAL) {
TRACE("prepare_cleanup_device");
@@ -2748,7 +2761,11 @@
}
if(mp->approximate_hess_kl){ cudaFree(mp->d_hess_kl_crust_mantle);}
}
- cudaFree(mp->d_rmass_crust_mantle);
+ if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
+ cudaFree(mp->d_rmassx_crust_mantle);
+ cudaFree(mp->d_rmassy_crust_mantle);
+ }
+ cudaFree(mp->d_rmassz_crust_mantle);
//------------------------------------------
// outer_core
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -29,65 +29,165 @@
nspec_actually,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore, &
gammaxstore,gammaystore,gammazstore, &
- iregion_code,nglob,rmass,rhostore,kappavstore, &
- nglob_oceans,rmass_ocean_load,NSPEC2D_TOP,ibelm_top,jacobian2D_top, &
- xstore,ystore,zstore,RHO_OCEANS)
+ iregion_code,rhostore,kappavstore, &
+ nglob_xy,nglob,prname, &
+ rmassx,rmassy,rmassz, &
+ nglob_oceans,rmass_ocean_load, &
+ xstore,ystore,zstore,RHO_OCEANS, &
+ NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax, &
+ rho_vp,rho_vs,nspec_stacey, &
+ jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
+ jacobian2D_bottom,jacobian2D_top)
-! creates rmass and rmass_ocean_load
+ ! creates rmassx, rmassy, rmassz and rmass_ocean_load
use meshfem3D_models_par
+ use meshfem3D_par,only: DT, NCHUNKS, ABSORBING_CONDITIONS, ichunk
implicit none
- integer myrank,nspec
+ integer :: myrank,nspec
+ integer :: idoubling(nspec)
+ integer :: ibool(NGLLX,NGLLY,NGLLZ,nspec)
+ integer :: nspec_actually
- integer idoubling(nspec)
-
- double precision wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ)
-
- integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
- integer nspec_actually
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_actually) :: &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore
- integer iregion_code
+ integer :: iregion_code,nglob_xy,nglob
- ! mass matrix
- integer nglob
- real(kind=CUSTOM_REAL), dimension(nglob) :: rmass
+ ! mass matrices
+ ! add C*deltat/2 contribution to the mass matrices on Stacey edges
+ real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rhostore,kappavstore
! ocean mass matrix
- integer nglob_oceans
+ integer :: nglob_oceans
real(kind=CUSTOM_REAL), dimension(nglob_oceans) :: rmass_ocean_load
- integer NSPEC2D_TOP
- integer, dimension(NSPEC2D_TOP) :: ibelm_top
+ ! arrays with the mesh in double precision
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ystore
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: zstore
+
+ double precision :: RHO_OCEANS
+
+ ! processor identification
+ character(len=150) prname
+
+ ! time scheme
+ real(kind=CUSTOM_REAL) :: deltat,deltatover2
+
+ ! Stacey conditions put back
+ integer :: NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
+ integer :: nspec_stacey
+
+ double precision, dimension(NGLLX) :: wxgll
+ double precision, dimension(NGLLY) :: wygll
+ double precision, dimension(NGLLZ) :: wzgll
+
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX) :: jacobian2D_xmin,jacobian2D_xmax
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX) :: jacobian2D_ymin,jacobian2D_ymax
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_BOTTOM) :: jacobian2D_bottom
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_TOP) :: jacobian2D_top
- ! arrays with the mesh in double precision
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX) :: normal_xmin,normal_xmax
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX) :: normal_ymin,normal_ymax
- double precision RHO_OCEANS
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_stacey) :: rho_vp,rho_vs
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+ real(kind=CUSTOM_REAL) :: tx,ty,tz,sn
+ real(kind=CUSTOM_REAL) :: nx,ny,nz,vn
+
+ integer, dimension(NSPEC2D_TOP) :: ibelm_top
+ integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
+ integer, dimension(NSPEC2DMAX_XMIN_XMAX) :: ibelm_xmin,ibelm_xmax
+ integer, dimension(NSPEC2DMAX_YMIN_YMAX) :: ibelm_ymin,ibelm_ymax
+
+ integer, dimension(2,NSPEC2DMAX_YMIN_YMAX) :: nimin,nimax,nkmin_eta
+ integer, dimension(2,NSPEC2DMAX_XMIN_XMAX) :: njmin,njmax,nkmin_xi
+
+ integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,ispec2D
+
! local parameters
- double precision weight
- double precision xval,yval,zval,rval,thetaval,phival
- double precision lat,lon,colat
- double precision elevation,height_oceans
+ double precision :: xval,yval,zval,rval,thetaval,phival,weight
+ double precision :: lat,lon,colat
+ double precision :: elevation,height_oceans
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
- integer :: ispec,i,j,k,iglobnum
+ integer :: ispec,i,j,k,iglob
integer :: ix_oceans,iy_oceans,iz_oceans,ispec_oceans,ispec2D_top_crust
+ ! initializes matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
- ! initializes
- rmass(:) = 0._CUSTOM_REAL
+ rmassx(:) = 0._CUSTOM_REAL
+ rmassy(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
+ ! use the non-dimensional time step to make the mass matrix correction
+ deltat = DT*dsqrt(PI*GRAV*RHOAV)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ do i=1,NGLLX
+ do j=1,NGLLY
+ wgllwgll_xy(i,j) = sngl(wxgll(i)*wygll(j))
+ enddo
+ enddo
+
+ do i=1,NGLLX
+ do k=1,NGLLZ
+ wgllwgll_xz(i,k) = sngl(wxgll(i)*wzgll(k))
+ enddo
+ enddo
+
+ do j=1,NGLLY
+ do k=1,NGLLZ
+ wgllwgll_yz(j,k) = sngl(wygll(j)*wzgll(k))
+ enddo
+ enddo
+
+ deltatover2 = sngl(0.5d0*deltat)
+
+ else ! double precision version
+ do i=1,NGLLX
+ do j=1,NGLLY
+ wgllwgll_xy(i,j) = wxgll(i)*wygll(j)
+ enddo
+ enddo
+
+ do i=1,NGLLX
+ do k=1,NGLLZ
+ wgllwgll_xz(i,k) = wxgll(i)*wzgll(k)
+ enddo
+ enddo
+
+ do j=1,NGLLY
+ do k=1,NGLLZ
+ wgllwgll_yz(j,k) = wygll(j)*wzgll(k)
+ enddo
+ enddo
+
+ deltatover2 = 0.5d0*deltat
+
+ endif
+
do ispec=1,nspec
! suppress fictitious elements in central cube
@@ -98,7 +198,7 @@
do i = 1,NGLLX
weight = wxgll(i)*wygll(j)*wzgll(k)
- iglobnum = ibool(i,j,k,ispec)
+ iglob = ibool(i,j,k,ispec)
! compute the jacobian
xixl = xixstore(i,j,k,ispec)
@@ -121,10 +221,10 @@
case( IREGION_CRUST_MANTLE, IREGION_INNER_CORE )
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglobnum) = rmass(iglobnum) + &
+ rmassz(iglob) = rmassz(iglob) + &
sngl(dble(rhostore(i,j,k,ispec)) * dble(jacobianl) * weight)
else
- rmass(iglobnum) = rmass(iglobnum) + rhostore(i,j,k,ispec) * jacobianl * weight
+ rmassz(iglob) = rmassz(iglob) + rhostore(i,j,k,ispec) * jacobianl * weight
endif
! fluid in outer core
@@ -134,10 +234,10 @@
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglobnum) = rmass(iglobnum) + &
+ rmassz(iglob) = rmassz(iglob) + &
sngl(dble(jacobianl) * weight * dble(rhostore(i,j,k,ispec)) / dble(kappavstore(i,j,k,ispec)))
else
- rmass(iglobnum) = rmass(iglobnum) + &
+ rmassz(iglob) = rmassz(iglob) + &
jacobianl * weight * rhostore(i,j,k,ispec) / kappavstore(i,j,k,ispec)
endif
@@ -168,7 +268,7 @@
do ix_oceans = 1,NGLLX
do iy_oceans = 1,NGLLY
- iglobnum=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
+ iglob=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
! if 3D Earth, compute local height of oceans
if(CASE_3D) then
@@ -213,9 +313,9 @@
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + sngl(weight)
+ rmass_ocean_load(iglob) = rmass_ocean_load(iglob) + sngl(weight)
else
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + weight
+ rmass_ocean_load(iglob) = rmass_ocean_load(iglob) + weight
endif
enddo
@@ -224,8 +324,346 @@
enddo
! add regular mass matrix to ocean load contribution
- rmass_ocean_load(:) = rmass_ocean_load(:) + rmass(:)
+ rmass_ocean_load(:) = rmass_ocean_load(:) + rmassz(:)
endif
+ ! add C*deltat/2 contribution to the mass matrices on Stacey edges
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+
+ ! read arrays for Stacey conditions
+ open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
+ status='old',form='unformatted',action='read')
+ read(27) nimin
+ read(27) nimax
+ read(27) njmin
+ read(27) njmax
+ read(27) nkmin_xi
+ read(27) nkmin_eta
+ close(27)
+
+ select case(iregion_code)
+
+ case(IREGION_CRUST_MANTLE)
+
+ rmassx(:) = rmassz(:)
+ rmassy(:) = rmassz(:)
+
+ ! xmin
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AC) then
+
+ do ispec2D=1,nspec2D_xmin
+
+ ispec=ibelm_xmin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_xmin(1,j,k,ispec2D)
+ ny = normal_xmin(2,j,k,ispec2D)
+ nz = normal_xmin(3,j,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AC
+
+ ! xmax
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AB) then
+
+ do ispec2D=1,nspec2D_xmax
+
+ ispec=ibelm_xmax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_xmax(1,j,k,ispec2D)
+ ny = normal_xmax(2,j,k,ispec2D)
+ nz = normal_xmax(3,j,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AB
+
+ ! ymin
+ do ispec2D=1,nspec2D_ymin
+
+ ispec=ibelm_ymin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_ymin(1,i,k,ispec2D)
+ ny = normal_ymin(2,i,k,ispec2D)
+ nz = normal_ymin(3,i,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! ymax
+ do ispec2D=1,nspec2D_ymax
+
+ ispec=ibelm_ymax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_ymax(1,i,k,ispec2D)
+ ny = normal_ymax(2,i,k,ispec2D)
+ nz = normal_ymax(3,i,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! check that mass matrix is positive
+ if(minval(rmassx(:)) <= 0.) call exit_MPI(myrank,'negative rmassx matrix term')
+ if(minval(rmassy(:)) <= 0.) call exit_MPI(myrank,'negative rmassy matrix term')
+
+ case(IREGION_OUTER_CORE)
+
+ ! xmin
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AC) then
+
+ do ispec2D=1,nspec2D_xmin
+
+ ispec=ibelm_xmin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AC
+
+ ! xmax
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AB) then
+
+ do ispec2D=1,nspec2D_xmax
+
+ ispec=ibelm_xmax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AB
+
+ ! ymin
+ do ispec2D=1,nspec2D_ymin
+
+ ispec=ibelm_ymin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! ymax
+ do ispec2D=1,nspec2D_ymax
+
+ ispec=ibelm_ymax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! bottom (zmin)
+ do ispec2D=1,NSPEC2D_BOTTOM
+
+ ispec=ibelm_bottom(ispec2D)
+
+ k=1
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ case( IREGION_INNER_CORE )
+
+ case default
+ call exit_MPI(myrank,'wrong region code')
+
+ end select
+
+ endif
+
+ ! check that mass matrix is positive
+ if(minval(rmassz(:)) <= 0.) call exit_MPI(myrank,'negative rmassz matrix term')
+
end subroutine create_mass_matrices
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -189,8 +189,9 @@
integer :: nglob
integer :: ieoff,ilocnum,ier
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
+ ! mass matrices
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
+ integer :: nglob_xy
! mass matrix and bathymetry for ocean load
integer nglob_oceans
@@ -720,7 +721,7 @@
if(NCHUNKS /= 6) &
call get_absorb(myrank,prname,iboun,nspec,nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
-
+
! create AVS or DX mesh data for the slices
if(SAVE_MESH_FILES) then
call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
@@ -789,9 +790,33 @@
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,&
xigll,yigll,zigll)
- ! allocates mass matrix in this slice (will be fully assembled in the solver)
- allocate(rmass(nglob),stat=ier)
+ ! allocates mass matrices in this slice (will be fully assembled in the solver)
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+ select case(iregion_code)
+ case( IREGION_CRUST_MANTLE )
+ nglob_xy = nglob
+ case( IREGION_INNER_CORE, IREGION_OUTER_CORE )
+ nglob_xy = 1
+ endselect
+ else
+ nglob_xy = 1
+ endif
+
+ allocate(rmassx(nglob_xy),stat=ier)
if(ier /= 0) stop 'error in allocate 21'
+ allocate(rmassy(nglob_xy),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+ allocate(rmassz(nglob),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+
! allocates ocean load mass matrix as well if oceans
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) then
nglob_oceans = nglob
@@ -802,14 +827,23 @@
allocate(rmass_ocean_load(nglob_oceans),stat=ier)
if(ier /= 0) stop 'error in allocate 22'
- ! creating mass matrix in this slice (will be fully assembled in the solver)
+ ! creating mass matrices in this slice (will be fully assembled in the solver)
call create_mass_matrices(myrank,nspec,idoubling,wxgll,wygll,wzgll,ibool, &
nspec_actually,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore, &
gammaxstore,gammaystore,gammazstore, &
- iregion_code,nglob,rmass,rhostore,kappavstore, &
- nglob_oceans,rmass_ocean_load,NSPEC2D_TOP,ibelm_top,jacobian2D_top, &
- xstore,ystore,zstore,RHO_OCEANS)
+ iregion_code,rhostore,kappavstore, &
+ nglob_xy,nglob,prname, &
+ rmassx,rmassy,rmassz, &
+ nglob_oceans,rmass_ocean_load, &
+ xstore,ystore,zstore,RHO_OCEANS, &
+ NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax, &
+ rho_vp,rho_vs,nspec_stacey, &
+ jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
+ jacobian2D_bottom,jacobian2D_top)
! user output
if(myrank == 0 ) write(IMAIN,*) ' saving binary files'
@@ -823,12 +857,13 @@
nspec_ani,c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,is_on_a_slice_edge,rmass,rmass_ocean_load,nglob_oceans, &
+ ibool,idoubling,is_on_a_slice_edge,nglob_xy,nglob, &
+ rmassx,rmassy,rmassz,rmass_ocean_load,nglob_oceans, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top, &
jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top,nspec,nglob, &
+ jacobian2D_bottom,jacobian2D_top,nspec, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,OCEANS, &
@@ -836,7 +871,8 @@
size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4),size(tau_e_store,5),&
ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
- deallocate(rmass,rmass_ocean_load)
+ deallocate(rmassx,rmassy,rmassz)
+ deallocate(rmass_ocean_load)
! boundary mesh
if (SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -33,18 +33,20 @@
nspec_ani,c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,is_on_a_slice_edge,rmass,rmass_ocean_load,npointot_oceans, &
+ ibool,idoubling,is_on_a_slice_edge,nglob_xy,nglob, &
+ rmassx,rmassy,rmassz,rmass_ocean_load,npointot_oceans, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top, &
jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top,nspec,nglob, &
+ jacobian2D_bottom,jacobian2D_top,nspec, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION,vx,vy,vz,vnspec, &
ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
+ use meshfem3D_par,only: NCHUNKS
implicit none
@@ -55,7 +57,7 @@
character(len=150) prname
integer iregion_code
- integer nspec,nglob,nspec_stacey
+ integer nspec,nglob_xy,nglob,nspec_stacey
integer npointot_oceans
! Stacey
@@ -92,8 +94,9 @@
! this for non blocking MPI
logical, dimension(nspec) :: is_on_a_slice_edge
- ! mass matrix
- real(kind=CUSTOM_REAL) rmass(nglob)
+ ! mass matrices
+ real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
! additional ocean load mass matrix
real(kind=CUSTOM_REAL) rmass_ocean_load(npointot_oceans)
@@ -238,8 +241,20 @@
endif
- ! mass matrix
- write(27) rmass
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) then
+ write(27) rmassx
+ write(27) rmassy
+ endif
+
+ write(27) rmassz
! additional ocean load mass matrix if oceans and if we are in the crust
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) write(27) rmass_ocean_load
@@ -252,10 +267,10 @@
! mesh arrays used in the solver to locate source and receivers
! and for anisotropy and gravity, save in single precision
-! use rmass for temporary storage to perform conversion, since already saved
+! use rmassz for temporary storage to perform conversion, since already saved
!--- x coordinate
- rmass(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -263,18 +278,18 @@
iglob = ibool(i,j,k,ispec)
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglob) = sngl(xstore(i,j,k,ispec))
+ rmassz(iglob) = sngl(xstore(i,j,k,ispec))
else
- rmass(iglob) = xstore(i,j,k,ispec)
+ rmassz(iglob) = xstore(i,j,k,ispec)
endif
enddo
enddo
enddo
enddo
- write(27) rmass
+ write(27) rmassz
!--- y coordinate
- rmass(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -282,18 +297,18 @@
iglob = ibool(i,j,k,ispec)
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglob) = sngl(ystore(i,j,k,ispec))
+ rmassz(iglob) = sngl(ystore(i,j,k,ispec))
else
- rmass(iglob) = ystore(i,j,k,ispec)
+ rmassz(iglob) = ystore(i,j,k,ispec)
endif
enddo
enddo
enddo
enddo
- write(27) rmass
+ write(27) rmassz
!--- z coordinate
- rmass(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -301,15 +316,15 @@
iglob = ibool(i,j,k,ispec)
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglob) = sngl(zstore(i,j,k,ispec))
+ rmassz(iglob) = sngl(zstore(i,j,k,ispec))
else
- rmass(iglob) = zstore(i,j,k,ispec)
+ rmassz(iglob) = zstore(i,j,k,ispec)
endif
enddo
enddo
enddo
enddo
- write(27) rmass
+ write(27) rmassz
write(27) ibool
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -438,22 +438,37 @@
!
subroutine compute_coupling_ocean(accel_crust_mantle,b_accel_crust_mantle, &
- rmass_crust_mantle,rmass_ocean_load,normal_top_crust_mantle, &
+ rmassx_crust_mantle, rmassy_crust_mantle, rmassz_crust_mantle, &
+ rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load, &
- SIMULATION_TYPE,nspec_top)
+ updated_dof_ocean_load,NGLOB_XY, &
+ SIMULATION_TYPE,nspec_top, &
+ ABSORBING_CONDITIONS)
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
+ integer :: NGLOB_XY
+
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
accel_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE_ADJOINT) :: &
b_accel_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
+
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_CM) :: normal_top_crust_mantle
@@ -461,12 +476,15 @@
integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
logical, dimension(NGLOB_CRUST_MANTLE_OCEANS) :: updated_dof_ocean_load
+ logical :: ABSORBING_CONDITIONS
integer SIMULATION_TYPE
integer nspec_top
! local parameters
real(kind=CUSTOM_REAL) :: force_normal_comp,b_force_normal_comp
+ real(kind=CUSTOM_REAL) :: additional_term_x,additional_term_y,additional_term_z
+ real(kind=CUSTOM_REAL) :: b_additional_term_x,b_additional_term_y,b_additional_term_z
real(kind=CUSTOM_REAL) :: additional_term,b_additional_term
real(kind=CUSTOM_REAL) :: nx,ny,nz
integer :: i,j,k,ispec,ispec2D,iglob
@@ -474,61 +492,127 @@
! initialize the updates
updated_dof_ocean_load(:) = .false.
- ! for surface elements exactly at the top of the crust (ocean bottom)
- do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
- ispec = ibelm_top_crust_mantle(ispec2D)
+ ! for surface elements exactly at the top of the crust (ocean bottom)
+ do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
- ! only for DOFs exactly at the top of the crust (ocean bottom)
- k = NGLLZ
+ ispec = ibelm_top_crust_mantle(ispec2D)
- do j = 1,NGLLY
- do i = 1,NGLLX
+ ! only for DOFs exactly at the top of the crust (ocean bottom)
+ k = NGLLZ
- ! get global point number
- iglob = ibool_crust_mantle(i,j,k,ispec)
+ do j = 1,NGLLY
+ do i = 1,NGLLX
- ! only update once
- if(.not. updated_dof_ocean_load(iglob)) then
+ ! get global point number
+ iglob = ibool_crust_mantle(i,j,k,ispec)
- ! get normal
- nx = normal_top_crust_mantle(1,i,j,ispec2D)
- ny = normal_top_crust_mantle(2,i,j,ispec2D)
- nz = normal_top_crust_mantle(3,i,j,ispec2D)
+ ! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
- ! make updated component of right-hand side
- ! we divide by rmass_crust_mantle() which is 1 / M
- ! we use the total force which includes the Coriolis term above
- force_normal_comp = (accel_crust_mantle(1,iglob)*nx + &
- accel_crust_mantle(2,iglob)*ny + &
- accel_crust_mantle(3,iglob)*nz) / rmass_crust_mantle(iglob)
+ ! get normal
+ nx = normal_top_crust_mantle(1,i,j,ispec2D)
+ ny = normal_top_crust_mantle(2,i,j,ispec2D)
+ nz = normal_top_crust_mantle(3,i,j,ispec2D)
- additional_term = (rmass_ocean_load(iglob) - rmass_crust_mantle(iglob)) * force_normal_comp
+ ! make updated component of right-hand side
+ ! we divide by rmass_crust_mantle() which is 1 / M
+ ! we use the total force which includes the Coriolis term above
+ force_normal_comp = accel_crust_mantle(1,iglob)*nx / rmassx_crust_mantle(iglob) + &
+ accel_crust_mantle(2,iglob)*ny / rmassy_crust_mantle(iglob) + &
+ accel_crust_mantle(3,iglob)*nz / rmassz_crust_mantle(iglob)
- accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term * nx
- accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term * ny
- accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term * nz
+ additional_term_x = (rmass_ocean_load(iglob) - rmassx_crust_mantle(iglob)) * force_normal_comp
+ additional_term_y = (rmass_ocean_load(iglob) - rmassy_crust_mantle(iglob)) * force_normal_comp
+ additional_term_z = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
- if (SIMULATION_TYPE == 3) then
- b_force_normal_comp = (b_accel_crust_mantle(1,iglob)*nx + &
- b_accel_crust_mantle(2,iglob)*ny + &
- b_accel_crust_mantle(3,iglob)*nz) / rmass_crust_mantle(iglob)
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term_x * nx
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term_y * ny
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term_z * nz
- b_additional_term = (rmass_ocean_load(iglob) - rmass_crust_mantle(iglob)) * b_force_normal_comp
+ if (SIMULATION_TYPE == 3) then
+ b_force_normal_comp = b_accel_crust_mantle(1,iglob)*nx / rmassx_crust_mantle(iglob) + &
+ b_accel_crust_mantle(2,iglob)*ny / rmassy_crust_mantle(iglob) + &
+ b_accel_crust_mantle(3,iglob)*nz / rmassz_crust_mantle(iglob)
- b_accel_crust_mantle(1,iglob) = b_accel_crust_mantle(1,iglob) + b_additional_term * nx
- b_accel_crust_mantle(2,iglob) = b_accel_crust_mantle(2,iglob) + b_additional_term * ny
- b_accel_crust_mantle(3,iglob) = b_accel_crust_mantle(3,iglob) + b_additional_term * nz
- endif
+ b_additional_term_x = (rmass_ocean_load(iglob) - rmassx_crust_mantle(iglob)) * b_force_normal_comp
+ b_additional_term_y = (rmass_ocean_load(iglob) - rmassy_crust_mantle(iglob)) * b_force_normal_comp
+ b_additional_term_z = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * b_force_normal_comp
- ! done with this point
- updated_dof_ocean_load(iglob) = .true.
+ b_accel_crust_mantle(1,iglob) = b_accel_crust_mantle(1,iglob) + b_additional_term_x * nx
+ b_accel_crust_mantle(2,iglob) = b_accel_crust_mantle(2,iglob) + b_additional_term_y * ny
+ b_accel_crust_mantle(3,iglob) = b_accel_crust_mantle(3,iglob) + b_additional_term_z * nz
+ endif
- endif
+ ! done with this point
+ updated_dof_ocean_load(iglob) = .true.
- enddo
- enddo
- enddo
+ endif
+ enddo
+ enddo
+ enddo
+
+ else
+
+ ! for surface elements exactly at the top of the crust (ocean bottom)
+ do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+
+ ispec = ibelm_top_crust_mantle(ispec2D)
+
+ ! only for DOFs exactly at the top of the crust (ocean bottom)
+ k = NGLLZ
+
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ ! get global point number
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+
+ ! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
+
+ ! get normal
+ nx = normal_top_crust_mantle(1,i,j,ispec2D)
+ ny = normal_top_crust_mantle(2,i,j,ispec2D)
+ nz = normal_top_crust_mantle(3,i,j,ispec2D)
+
+ ! make updated component of right-hand side
+ ! we divide by rmass_crust_mantle() which is 1 / M
+ ! we use the total force which includes the Coriolis term above
+ force_normal_comp = (accel_crust_mantle(1,iglob)*nx + &
+ accel_crust_mantle(2,iglob)*ny + &
+ accel_crust_mantle(3,iglob)*nz) / rmassz_crust_mantle(iglob)
+
+ additional_term = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
+
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term * nx
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term * ny
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term * nz
+
+ if (SIMULATION_TYPE == 3) then
+ b_force_normal_comp = (b_accel_crust_mantle(1,iglob)*nx + &
+ b_accel_crust_mantle(2,iglob)*ny + &
+ b_accel_crust_mantle(3,iglob)*nz) / rmassz_crust_mantle(iglob)
+
+ b_additional_term = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * b_force_normal_comp
+
+ b_accel_crust_mantle(1,iglob) = b_accel_crust_mantle(1,iglob) + b_additional_term * nx
+ b_accel_crust_mantle(2,iglob) = b_accel_crust_mantle(2,iglob) + b_additional_term * ny
+ b_accel_crust_mantle(3,iglob) = b_accel_crust_mantle(3,iglob) + b_additional_term * nz
+ endif
+
+ ! done with this point
+ updated_dof_ocean_load(iglob) = .true.
+
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ endif
+
end subroutine compute_coupling_ocean
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -45,7 +45,11 @@
epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ tempx1_att,tempx2_att,tempx3_att, &
+ tempy1_att,tempy2_att,tempy3_att, &
+ tempz1_att,tempz2_att,tempz3_att, &
+ epsilondev_loc,rho_s_H)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -99,6 +103,9 @@
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
@@ -111,6 +118,10 @@
real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
real(kind=CUSTOM_REAL) templ
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
@@ -170,20 +181,57 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- epsilondev_loc(1,i,j,k) = duxdxl - templ
- epsilondev_loc(2,i,j,k) = duydyl - templ
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! 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
+
+ ! compute deviatoric strain
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
+ endif
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilondev_loc(1,i,j,k) = duxdxl_att - templ
+ epsilondev_loc(2,i,j,k) = duydyl_att - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ !$OMP CRITICAL
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ !$OMP END CRITICAL
+ else
+ ispec_strain = ispec
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ endif
+ epsilondev_loc(1,i,j,k) = duxdxl - templ
+ epsilondev_loc(2,i,j,k) = duydyl - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
! precompute terms for attenuation if needed
@@ -377,7 +425,11 @@
epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ tempx1_att,tempx2_att,tempx3_att, &
+ tempy1_att,tempy2_att,tempy3_att, &
+ tempz1_att,tempz2_att,tempz3_att, &
+ epsilondev_loc,rho_s_H)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -435,6 +487,10 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
@@ -460,6 +516,10 @@
real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
real(kind=CUSTOM_REAL) templ
real(kind=CUSTOM_REAL) templ1,templ1_cos,templ2,templ2_cos,templ3,templ3_two,templ3_cos
real(kind=CUSTOM_REAL) kappavl,kappahl,muvl,muhl
@@ -520,20 +580,57 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- epsilondev_loc(1,i,j,k) = duxdxl - templ
- epsilondev_loc(2,i,j,k) = duydyl - templ
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! 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
+
+ ! compute deviatoric strain
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
+ endif
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilondev_loc(1,i,j,k) = duxdxl_att - templ
+ epsilondev_loc(2,i,j,k) = duydyl_att - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ !$OMP CRITICAL
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ !$OMP END CRITICAL
+ else
+ ispec_strain = ispec
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ endif
+ epsilondev_loc(1,i,j,k) = duxdxl - templ
+ epsilondev_loc(2,i,j,k) = duydyl - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
! precompute terms for attenuation if needed
@@ -919,7 +1016,11 @@
epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ tempx1_att,tempx2_att,tempx3_att, &
+ tempy1_att,tempy2_att,tempy3_att, &
+ tempz1_att,tempz2_att,tempz3_att, &
+ epsilondev_loc,rho_s_H)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -975,6 +1076,10 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
+
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
@@ -990,6 +1095,10 @@
real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) duxdxl_att,duxdyl_att,duxdzl_att,duydxl_att
+ real(kind=CUSTOM_REAL) duydyl_att,duydzl_att,duzdxl_att,duzdyl_att,duzdzl_att
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl_att,duzdxl_plus_duxdzl_att,duzdyl_plus_duydzl_att
+
real(kind=CUSTOM_REAL) templ
! for gravity
@@ -1048,20 +1157,57 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- epsilondev_loc(1,i,j,k) = duxdxl - templ
- epsilondev_loc(2,i,j,k) = duydyl - templ
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ if ( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
+ ! temporary variables used for fixing attenuation in a consistent way
+ duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ duxdzl_att = xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
+
+ duydxl_att = xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ duydyl_att = xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ duydzl_att = xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
+
+ duzdxl_att = xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ duzdyl_att = xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ duzdzl_att = xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+
+ ! 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
+
+ ! compute deviatoric strain
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
+ endif
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilondev_loc(1,i,j,k) = duxdxl_att - templ
+ epsilondev_loc(2,i,j,k) = duydyl_att - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ else
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ !$OMP CRITICAL
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ !$OMP END CRITICAL
+ else
+ ispec_strain = ispec
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ endif
+ epsilondev_loc(1,i,j,k) = duxdxl - templ
+ epsilondev_loc(2,i,j,k) = duydyl - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
endif
! precompute terms for attenuation if needed
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic.F90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -273,14 +273,16 @@
! on CPU
call compute_forces_ac_update_veloc(NGLOB_OUTER_CORE,veloc_outer_core,accel_outer_core, &
deltatover2,rmass_outer_core)
- ! adjoint / kernel runs
+
+ ! adjoint / kernel runs
if (SIMULATION_TYPE == 3) &
call compute_forces_ac_update_veloc(NGLOB_OUTER_CORE_ADJOINT,b_veloc_outer_core,b_accel_outer_core, &
b_deltatover2,rmass_outer_core)
+
else
! on GPU
call kernel_3_outer_core_cuda(Mesh_pointer, &
- deltatover2,SIMULATION_TYPE,b_deltatover2)
+ deltatover2,SIMULATION_TYPE,b_deltatover2)
endif
end subroutine compute_forces_acoustic
@@ -288,7 +290,7 @@
!=====================================================================
subroutine compute_forces_ac_update_veloc(NGLOB,veloc_outer_core,accel_outer_core, &
- deltatover2,rmass_outer_core)
+ deltatover2,rmass_outer_core)
use constants_solver,only: CUSTOM_REAL
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -231,7 +231,7 @@
tempz3l = tempz3l + displ_crust_mantle(3,iglob)*hp3
enddo
- if(ATTENUATION_VAL) then
+ if( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
! temporary variables used for fixing attenuation in a consistent way
tempx1l_att = 0._CUSTOM_REAL
@@ -336,7 +336,7 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- if(ATTENUATION_VAL) then
+ if( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
! 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
@@ -356,19 +356,17 @@
duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
- epsilondev_loc(1,i,j,k) = duxdxl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(2,i,j,k) = duydyl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
endif
+ epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilondev_loc(1,i,j,k) = duxdxl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(2,i,j,k) = duydyl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
else
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -256,7 +256,7 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
if(ATTENUATION_NEW_VAL) then
! takes new routines
! use first order Taylor expansion of displacement for local storage of stresses
@@ -270,72 +270,72 @@
! since we know that NGLLX = 5, this should help pipelining
iglobv5(:) = ibool(:,j,k,ispec)
- dummyx_loc_att(1,j,k) = displ_crust_mantle(1,iglobv5(1)) + deltat*veloc_crust_mantle(1,iglobv5(1))
- dummyy_loc_att(1,j,k) = displ_crust_mantle(2,iglobv5(1)) + deltat*veloc_crust_mantle(2,iglobv5(1))
- dummyz_loc_att(1,j,k) = displ_crust_mantle(3,iglobv5(1)) + deltat*veloc_crust_mantle(3,iglobv5(1))
+ dummyx_loc_att(1,j,k) = deltat*veloc_crust_mantle(1,iglobv5(1))
+ dummyy_loc_att(1,j,k) = deltat*veloc_crust_mantle(2,iglobv5(1))
+ dummyz_loc_att(1,j,k) = deltat*veloc_crust_mantle(3,iglobv5(1))
- dummyx_loc_att(2,j,k) = displ_crust_mantle(1,iglobv5(2)) + deltat*veloc_crust_mantle(1,iglobv5(2))
- dummyy_loc_att(2,j,k) = displ_crust_mantle(2,iglobv5(2)) + deltat*veloc_crust_mantle(2,iglobv5(2))
- dummyz_loc_att(2,j,k) = displ_crust_mantle(3,iglobv5(2)) + deltat*veloc_crust_mantle(3,iglobv5(2))
+ dummyx_loc_att(2,j,k) = deltat*veloc_crust_mantle(1,iglobv5(2))
+ dummyy_loc_att(2,j,k) = deltat*veloc_crust_mantle(2,iglobv5(2))
+ dummyz_loc_att(2,j,k) = deltat*veloc_crust_mantle(3,iglobv5(2))
- dummyx_loc_att(3,j,k) = displ_crust_mantle(1,iglobv5(3)) + deltat*veloc_crust_mantle(1,iglobv5(3))
- dummyy_loc_att(3,j,k) = displ_crust_mantle(2,iglobv5(3)) + deltat*veloc_crust_mantle(2,iglobv5(3))
- dummyz_loc_att(3,j,k) = displ_crust_mantle(3,iglobv5(3)) + deltat*veloc_crust_mantle(3,iglobv5(3))
+ dummyx_loc_att(3,j,k) = deltat*veloc_crust_mantle(1,iglobv5(3))
+ dummyy_loc_att(3,j,k) = deltat*veloc_crust_mantle(2,iglobv5(3))
+ dummyz_loc_att(3,j,k) = deltat*veloc_crust_mantle(3,iglobv5(3))
- dummyx_loc_att(4,j,k) = displ_crust_mantle(1,iglobv5(4)) + deltat*veloc_crust_mantle(1,iglobv5(4))
- dummyy_loc_att(4,j,k) = displ_crust_mantle(2,iglobv5(4)) + deltat*veloc_crust_mantle(2,iglobv5(4))
- dummyz_loc_att(4,j,k) = displ_crust_mantle(3,iglobv5(4)) + deltat*veloc_crust_mantle(3,iglobv5(4))
+ dummyx_loc_att(4,j,k) = deltat*veloc_crust_mantle(1,iglobv5(4))
+ dummyy_loc_att(4,j,k) = deltat*veloc_crust_mantle(2,iglobv5(4))
+ dummyz_loc_att(4,j,k) = deltat*veloc_crust_mantle(3,iglobv5(4))
- dummyx_loc_att(5,j,k) = displ_crust_mantle(1,iglobv5(5)) + deltat*veloc_crust_mantle(1,iglobv5(5))
- dummyy_loc_att(5,j,k) = displ_crust_mantle(2,iglobv5(5)) + deltat*veloc_crust_mantle(2,iglobv5(5))
- dummyz_loc_att(5,j,k) = displ_crust_mantle(3,iglobv5(5)) + deltat*veloc_crust_mantle(3,iglobv5(5))
+ dummyx_loc_att(5,j,k) = deltat*veloc_crust_mantle(1,iglobv5(5))
+ dummyy_loc_att(5,j,k) = deltat*veloc_crust_mantle(2,iglobv5(5))
+ dummyz_loc_att(5,j,k) = deltat*veloc_crust_mantle(3,iglobv5(5))
#else
! way 1:
do i=1,NGLLX
iglob1 = ibool(i,j,k,ispec)
- dummyx_loc_att(i,j,k) = displ_crust_mantle(1,iglob1) + deltat*veloc_crust_mantle(1,iglob1)
- dummyy_loc_att(i,j,k) = displ_crust_mantle(2,iglob1) + deltat*veloc_crust_mantle(2,iglob1)
- dummyz_loc_att(i,j,k) = displ_crust_mantle(3,iglob1) + deltat*veloc_crust_mantle(3,iglob1)
+ dummyx_loc_att(i,j,k) = deltat*veloc_crust_mantle(1,iglob1)
+ dummyy_loc_att(i,j,k) = deltat*veloc_crust_mantle(2,iglob1)
+ dummyz_loc_att(i,j,k) = deltat*veloc_crust_mantle(3,iglob1)
enddo
#endif
enddo
enddo
- endif
- else
- ! takes old routines
- do k=1,NGLLZ
- do j=1,NGLLY
+ else
+ ! takes old routines
+ do k=1,NGLLZ
+ do j=1,NGLLY
#ifdef _HANDOPT
- dummyx_loc_att(1,j,k) = dummyx_loc(1,j,k)
- dummyy_loc_att(1,j,k) = dummyx_loc(1,j,k)
- dummyz_loc_att(1,j,k) = dummyx_loc(1,j,k)
-
- dummyx_loc_att(2,j,k) = dummyx_loc(2,j,k)
- dummyy_loc_att(2,j,k) = dummyx_loc(2,j,k)
- dummyz_loc_att(2,j,k) = dummyx_loc(2,j,k)
-
- dummyx_loc_att(3,j,k) = dummyx_loc(3,j,k)
- dummyy_loc_att(3,j,k) = dummyx_loc(3,j,k)
- dummyz_loc_att(3,j,k) = dummyx_loc(3,j,k)
-
- dummyx_loc_att(4,j,k) = dummyx_loc(4,j,k)
- dummyy_loc_att(4,j,k) = dummyx_loc(4,j,k)
- dummyz_loc_att(4,j,k) = dummyx_loc(4,j,k)
-
- dummyx_loc_att(5,j,k) = dummyx_loc(5,j,k)
- dummyy_loc_att(5,j,k) = dummyx_loc(5,j,k)
- dummyz_loc_att(5,j,k) = dummyx_loc(5,j,k)
+ dummyx_loc_att(1,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(1,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(1,j,k) = 0._CUSTOM_REAL
+
+ dummyx_loc_att(2,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(2,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(2,j,k) = 0._CUSTOM_REAL
+
+ dummyx_loc_att(3,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(3,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(3,j,k) = 0._CUSTOM_REAL
+
+ dummyx_loc_att(4,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(4,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(4,j,k) = 0._CUSTOM_REAL
+
+ dummyx_loc_att(5,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(5,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(5,j,k) = 0._CUSTOM_REAL
#else
- do i=1,NGLLX
- dummyx_loc_att(i,j,k) = dummyx_loc(i,j,k)
- dummyy_loc_att(i,j,k) = dummyy_loc(i,j,k)
- dummyz_loc_att(i,j,k) = dummyz_loc(i,j,k)
- enddo
+ do i=1,NGLLX
+ dummyx_loc_att(i,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(i,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(i,j,k) = 0._CUSTOM_REAL
+ enddo
#endif
+ enddo
enddo
- enddo
+ endif
endif
do j=1,m2
@@ -360,29 +360,36 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
do j=1,m2
do i=1,m1
- C1_m1_m2_5points_att(i,j) = hprime_xx(i,1)*B1_m1_m2_5points_att(1,j) + &
+ C1_m1_m2_5points_att(i,j) = C1_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_5points_att(1,j) + &
hprime_xx(i,2)*B1_m1_m2_5points_att(2,j) + &
hprime_xx(i,3)*B1_m1_m2_5points_att(3,j) + &
hprime_xx(i,4)*B1_m1_m2_5points_att(4,j) + &
hprime_xx(i,5)*B1_m1_m2_5points_att(5,j)
- C2_m1_m2_5points_att(i,j) = hprime_xx(i,1)*B2_m1_m2_5points_att(1,j) + &
+ C2_m1_m2_5points_att(i,j) = C2_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_5points_att(1,j) + &
hprime_xx(i,2)*B2_m1_m2_5points_att(2,j) + &
hprime_xx(i,3)*B2_m1_m2_5points_att(3,j) + &
hprime_xx(i,4)*B2_m1_m2_5points_att(4,j) + &
hprime_xx(i,5)*B2_m1_m2_5points_att(5,j)
- C3_m1_m2_5points_att(i,j) = hprime_xx(i,1)*B3_m1_m2_5points_att(1,j) + &
+ C3_m1_m2_5points_att(i,j) = C3_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_5points_att(1,j) + &
hprime_xx(i,2)*B3_m1_m2_5points_att(2,j) + &
hprime_xx(i,3)*B3_m1_m2_5points_att(3,j) + &
hprime_xx(i,4)*B3_m1_m2_5points_att(4,j) + &
hprime_xx(i,5)*B3_m1_m2_5points_att(5,j)
enddo
enddo
+ else
+ tempx1_att(:,:,:) = 0._CUSTOM_REAL
+ tempy1_att(:,:,:) = 0._CUSTOM_REAL
+ tempz1_att(:,:,:) = 0._CUSTOM_REAL
endif
do j=1,m1
@@ -410,25 +417,28 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
do j=1,m1
do i=1,m1
! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
do k = 1,NGLLX
- tempx2_att(i,j,k) = dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
dummyx_loc_att(i,5,k)*hprime_xxT(5,j)
- tempy2_att(i,j,k) = dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
dummyy_loc_att(i,5,k)*hprime_xxT(5,j)
- tempz2_att(i,j,k) = dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
@@ -436,6 +446,10 @@
enddo
enddo
enddo
+ else
+ tempx2_att(:,:,:) = 0._CUSTOM_REAL
+ tempy2_att(:,:,:) = 0._CUSTOM_REAL
+ tempz2_att(:,:,:) = 0._CUSTOM_REAL
endif
do j=1,m1
@@ -460,137 +474,103 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
do j=1,m1
do i=1,m2
- C1_mxm_m2_m1_5points_att(i,j) = A1_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ C1_mxm_m2_m1_5points_att(i,j) = C1_mxm_m2_m1_5points(i,j) + &
+ A1_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
A1_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
A1_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
A1_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
A1_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
- C2_mxm_m2_m1_5points_att(i,j) = A2_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ C2_mxm_m2_m1_5points_att(i,j) = C2_mxm_m2_m1_5points(i,j) + &
+ A2_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
A2_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
A2_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
A2_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
A2_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
- C3_mxm_m2_m1_5points_att(i,j) = A3_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ C3_mxm_m2_m1_5points_att(i,j) = C3_mxm_m2_m1_5points(i,j) + &
+ A3_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
A3_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
A3_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
A3_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
A3_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
enddo
enddo
+ else
+ tempx3_att(:,:,:) = 0._CUSTOM_REAL
+ tempy3_att(:,:,:) = 0._CUSTOM_REAL
+ tempz3_att(:,:,:) = 0._CUSTOM_REAL
endif
!
! compute either isotropic, transverse isotropic or anisotropic elements
!
if(ANISOTROPIC_3D_MANTLE_VAL) then
- ! anisotropic element
+ ! anisotropic element
- if(ATTENUATION_VAL) then
- call compute_element_aniso(ispec, &
+ call compute_element_aniso(ispec, &
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vx,vy,vz,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ tempx1_att,tempx2_att,tempx3_att, &
+ tempy1_att,tempy2_att,tempy3_att, &
+ tempz1_att,tempz2_att,tempz3_att, &
+ epsilondev_loc,rho_s_H)
+ else
+
+ if( .not. ispec_is_tiso(ispec) ) then
+ ! isotropic element
+
+ call compute_element_iso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
wgll_cube, &
- c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
- c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
- c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ kappavstore,muvstore, &
ibool, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
tempx1_att,tempx2_att,tempx3_att, &
tempy1_att,tempy2_att,tempy3_att, &
tempz1_att,tempz2_att,tempz3_att, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ epsilondev_loc,rho_s_H)
else
- call compute_element_aniso(ispec, &
+ ! transverse isotropic element
+
+ call compute_element_tiso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
wgll_cube, &
- c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
- c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
- c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
ibool, &
R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- endif
- else
- if( .not. ispec_is_tiso(ispec) ) then
- ! isotropic element
-
- if(ATTENUATION_VAL) then
- call compute_element_iso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- kappavstore,muvstore, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1_att,tempx2_att,tempx3_att, &
- tempy1_att,tempy2_att,tempy3_att, &
- tempz1_att,tempz2_att,tempz3_att, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- else
- call compute_element_iso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- kappavstore,muvstore, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- endif
-
- else
- ! transverse isotropic element
-
- if(ATTENUATION_VAL) then
- call compute_element_tiso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1_att,tempx2_att,tempx3_att, &
- tempy1_att,tempy2_att,tempy3_att, &
- tempz1_att,tempz2_att,tempz3_att, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- else
- call compute_element_tiso(ispec, &
- minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- wgll_cube, &
- kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
- ibool, &
- R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- endif
- endif ! .not. ispec_is_tiso
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ tempx1_att,tempx2_att,tempx3_att, &
+ tempy1_att,tempy2_att,tempy3_att, &
+ tempz1_att,tempz2_att,tempz3_att, &
+ epsilondev_loc,rho_s_H)
+ endif ! .not. ispec_is_tiso
endif
! subroutines adapted from Deville, Fischer and Mund, High-order methods
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -42,7 +42,7 @@
! non blocking MPI
! iphase: iphase = 1 is for computing outer elements in the crust_mantle and inner_core regions,
! iphase = 2 is for computing inner elements (former icall parameter)
- integer :: iphase
+ integer :: iphase,NGLOB_XY
logical :: phase_is_inner
! ****************************************************
@@ -69,152 +69,152 @@
if( .NOT. GPU_MODE ) then
- ! on CPU
+ ! on CPU
- ! compute internal forces in the solid regions
- ! note: for anisotropy and gravity, x y and z contain r theta and phi
- if( USE_DEVILLE_PRODUCTS_VAL ) then
- ! uses Deville (2002) optimizations
- ! crust/mantle region
- call compute_forces_crust_mantle_Dev( NSPEC_CRUST_MANTLE_STR_OR_ATT,NGLOB_CRUST_MANTLE, &
- NSPEC_CRUST_MANTLE_ATTENUAT, &
- deltat, &
- displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle, &
- phase_is_inner, &
- R_xx_crust_mantle,R_yy_crust_mantle,R_xy_crust_mantle, &
- R_xz_crust_mantle,R_yz_crust_mantle, &
- epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
- epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
- eps_trace_over_3_crust_mantle, &
- alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
- ! inner core region
- call compute_forces_inner_core_Dev( NSPEC_INNER_CORE_STR_OR_ATT,NGLOB_INNER_CORE, &
- NSPEC_INNER_CORE_ATTENUATION, &
- deltat, &
- displ_inner_core,veloc_inner_core,accel_inner_core, &
- phase_is_inner, &
- R_xx_inner_core,R_yy_inner_core,R_xy_inner_core,R_xz_inner_core,R_yz_inner_core, &
- epsilondev_xx_inner_core,epsilondev_yy_inner_core,epsilondev_xy_inner_core, &
- epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
- eps_trace_over_3_inner_core,&
- alphaval,betaval,gammaval, &
- factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
-
- else
- ! no Deville optimization
- ! crust/mantle region
- call compute_forces_crust_mantle( NSPEC_CRUST_MANTLE_STR_OR_ATT,NGLOB_CRUST_MANTLE, &
- NSPEC_CRUST_MANTLE_ATTENUAT, &
- deltat, &
- displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle, &
- phase_is_inner, &
- R_xx_crust_mantle,R_yy_crust_mantle,R_xy_crust_mantle, &
- R_xz_crust_mantle,R_yz_crust_mantle, &
- epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
- epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
- eps_trace_over_3_crust_mantle, &
- alphaval,betaval,gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
- ! inner core region
- call compute_forces_inner_core( NSPEC_INNER_CORE_STR_OR_ATT,NGLOB_INNER_CORE, &
- NSPEC_INNER_CORE_ATTENUATION, &
- deltat, &
- displ_inner_core,veloc_inner_core,accel_inner_core, &
- phase_is_inner, &
- R_xx_inner_core,R_yy_inner_core,R_xy_inner_core,R_xz_inner_core,R_yz_inner_core, &
- epsilondev_xx_inner_core,epsilondev_yy_inner_core,epsilondev_xy_inner_core, &
- epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
- eps_trace_over_3_inner_core,&
- alphaval,betaval,gammaval, &
- factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
-
- endif
-
- ! adjoint / kernel runs
- if (SIMULATION_TYPE == 3 ) then
- if( USE_DEVILLE_PRODUCTS_VAL ) then
+ ! compute internal forces in the solid regions
+ ! note: for anisotropy and gravity, x y and z contain r theta and phi
+ if( USE_DEVILLE_PRODUCTS_VAL ) then
! uses Deville (2002) optimizations
! crust/mantle region
- call compute_forces_crust_mantle_Dev( NSPEC_CRUST_MANTLE_ADJOINT,NGLOB_CRUST_MANTLE_ADJOINT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT, &
- deltat, &
- b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
- phase_is_inner, &
- b_R_xx_crust_mantle,b_R_yy_crust_mantle,b_R_xy_crust_mantle, &
- b_R_xz_crust_mantle,b_R_yz_crust_mantle, &
- b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle,&
- b_epsilondev_xy_crust_mantle, &
- b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle, &
- b_eps_trace_over_3_crust_mantle, &
- b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
+ call compute_forces_crust_mantle_Dev( NSPEC_CRUST_MANTLE_STR_OR_ATT,NGLOB_CRUST_MANTLE, &
+ NSPEC_CRUST_MANTLE_ATTENUAT, &
+ deltat, &
+ displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle, &
+ phase_is_inner, &
+ R_xx_crust_mantle,R_yy_crust_mantle,R_xy_crust_mantle, &
+ R_xz_crust_mantle,R_yz_crust_mantle, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
+ eps_trace_over_3_crust_mantle, &
+ alphaval,betaval,gammaval,factor_common_crust_mantle, &
+ size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
! inner core region
- call compute_forces_inner_core_Dev( NSPEC_INNER_CORE_ADJOINT,NGLOB_INNER_CORE_ADJOINT, &
- NSPEC_INNER_CORE_STR_AND_ATT, &
- deltat, &
- b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
- phase_is_inner, &
- b_R_xx_inner_core,b_R_yy_inner_core,b_R_xy_inner_core, &
- b_R_xz_inner_core,b_R_yz_inner_core, &
- b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core,b_epsilondev_xy_inner_core, &
- b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core, &
- b_eps_trace_over_3_inner_core,&
- b_alphaval,b_betaval,b_gammaval, &
- factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
+ call compute_forces_inner_core_Dev( NSPEC_INNER_CORE_STR_OR_ATT,NGLOB_INNER_CORE, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ deltat, &
+ displ_inner_core,veloc_inner_core,accel_inner_core, &
+ phase_is_inner, &
+ R_xx_inner_core,R_yy_inner_core,R_xy_inner_core,R_xz_inner_core,R_yz_inner_core, &
+ epsilondev_xx_inner_core,epsilondev_yy_inner_core,epsilondev_xy_inner_core, &
+ epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
+ eps_trace_over_3_inner_core,&
+ alphaval,betaval,gammaval, &
+ factor_common_inner_core, &
+ size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
- else
+ else
! no Deville optimization
! crust/mantle region
- call compute_forces_crust_mantle( NSPEC_CRUST_MANTLE_ADJOINT,NGLOB_CRUST_MANTLE_ADJOINT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT, &
- deltat, &
- b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
- phase_is_inner, &
- b_R_xx_crust_mantle,b_R_yy_crust_mantle,b_R_xy_crust_mantle, &
- b_R_xz_crust_mantle,b_R_yz_crust_mantle, &
- b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle,&
- b_epsilondev_xy_crust_mantle, &
- b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle, &
- b_eps_trace_over_3_crust_mantle, &
- b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
- size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
-
+ call compute_forces_crust_mantle( NSPEC_CRUST_MANTLE_STR_OR_ATT,NGLOB_CRUST_MANTLE, &
+ NSPEC_CRUST_MANTLE_ATTENUAT, &
+ deltat, &
+ displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle, &
+ phase_is_inner, &
+ R_xx_crust_mantle,R_yy_crust_mantle,R_xy_crust_mantle, &
+ R_xz_crust_mantle,R_yz_crust_mantle, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
+ eps_trace_over_3_crust_mantle, &
+ alphaval,betaval,gammaval,factor_common_crust_mantle, &
+ size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
! inner core region
- call compute_forces_inner_core( NSPEC_INNER_CORE_ADJOINT,NGLOB_INNER_CORE_ADJOINT, &
- NSPEC_INNER_CORE_STR_AND_ATT, &
- deltat, &
- b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
- phase_is_inner, &
- b_R_xx_inner_core,b_R_yy_inner_core,b_R_xy_inner_core, &
- b_R_xz_inner_core,b_R_yz_inner_core, &
- b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core,b_epsilondev_xy_inner_core, &
- b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core, &
- b_eps_trace_over_3_inner_core,&
- b_alphaval,b_betaval,b_gammaval, &
- factor_common_inner_core, &
- size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
- endif
- endif !SIMULATION_TYPE == 3
+ call compute_forces_inner_core( NSPEC_INNER_CORE_STR_OR_ATT,NGLOB_INNER_CORE, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ deltat, &
+ displ_inner_core,veloc_inner_core,accel_inner_core, &
+ phase_is_inner, &
+ R_xx_inner_core,R_yy_inner_core,R_xy_inner_core,R_xz_inner_core,R_yz_inner_core, &
+ epsilondev_xx_inner_core,epsilondev_yy_inner_core,epsilondev_xy_inner_core, &
+ epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
+ eps_trace_over_3_inner_core,&
+ alphaval,betaval,gammaval, &
+ factor_common_inner_core, &
+ size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
+ endif
+
+ ! adjoint / kernel runs
+ if (SIMULATION_TYPE == 3 ) then
+ if( USE_DEVILLE_PRODUCTS_VAL ) then
+ ! uses Deville (2002) optimizations
+ ! crust/mantle region
+ call compute_forces_crust_mantle_Dev( NSPEC_CRUST_MANTLE_ADJOINT,NGLOB_CRUST_MANTLE_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT, &
+ deltat, &
+ b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ phase_is_inner, &
+ b_R_xx_crust_mantle,b_R_yy_crust_mantle,b_R_xy_crust_mantle, &
+ b_R_xz_crust_mantle,b_R_yz_crust_mantle, &
+ b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle,&
+ b_epsilondev_xy_crust_mantle, &
+ b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle, &
+ b_eps_trace_over_3_crust_mantle, &
+ b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
+ size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
+ ! inner core region
+ call compute_forces_inner_core_Dev( NSPEC_INNER_CORE_ADJOINT,NGLOB_INNER_CORE_ADJOINT, &
+ NSPEC_INNER_CORE_STR_AND_ATT, &
+ deltat, &
+ b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
+ phase_is_inner, &
+ b_R_xx_inner_core,b_R_yy_inner_core,b_R_xy_inner_core, &
+ b_R_xz_inner_core,b_R_yz_inner_core, &
+ b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core,b_epsilondev_xy_inner_core, &
+ b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core, &
+ b_eps_trace_over_3_inner_core,&
+ b_alphaval,b_betaval,b_gammaval, &
+ factor_common_inner_core, &
+ size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
+
+ else
+ ! no Deville optimization
+ ! crust/mantle region
+ call compute_forces_crust_mantle( NSPEC_CRUST_MANTLE_ADJOINT,NGLOB_CRUST_MANTLE_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT, &
+ deltat, &
+ b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ phase_is_inner, &
+ b_R_xx_crust_mantle,b_R_yy_crust_mantle,b_R_xy_crust_mantle, &
+ b_R_xz_crust_mantle,b_R_yz_crust_mantle, &
+ b_epsilondev_xx_crust_mantle,b_epsilondev_yy_crust_mantle,&
+ b_epsilondev_xy_crust_mantle, &
+ b_epsilondev_xz_crust_mantle,b_epsilondev_yz_crust_mantle, &
+ b_eps_trace_over_3_crust_mantle, &
+ b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
+ size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
+
+ ! inner core region
+ call compute_forces_inner_core( NSPEC_INNER_CORE_ADJOINT,NGLOB_INNER_CORE_ADJOINT, &
+ NSPEC_INNER_CORE_STR_AND_ATT, &
+ deltat, &
+ b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
+ phase_is_inner, &
+ b_R_xx_inner_core,b_R_yy_inner_core,b_R_xy_inner_core, &
+ b_R_xz_inner_core,b_R_yz_inner_core, &
+ b_epsilondev_xx_inner_core,b_epsilondev_yy_inner_core,b_epsilondev_xy_inner_core, &
+ b_epsilondev_xz_inner_core,b_epsilondev_yz_inner_core, &
+ b_eps_trace_over_3_inner_core,&
+ b_alphaval,b_betaval,b_gammaval, &
+ factor_common_inner_core, &
+ size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5) )
+ endif
+ endif !SIMULATION_TYPE == 3
+
else
- ! on GPU
- ! contains both forward SIM_TYPE==1 and backward SIM_TYPE==3 simulations
- ! for crust/mantle
- call compute_forces_crust_mantle_cuda(Mesh_pointer,deltat,iphase)
- ! for inner core
- call compute_forces_inner_core_cuda(Mesh_pointer,deltat,iphase)
+ ! on GPU
+ ! contains both forward SIM_TYPE==1 and backward SIM_TYPE==3 simulations
+ ! for crust/mantle
+ call compute_forces_crust_mantle_cuda(Mesh_pointer,deltat,iphase)
+ ! for inner core
+ call compute_forces_inner_core_cuda(Mesh_pointer,deltat,iphase)
endif ! GPU_MODE
@@ -491,19 +491,29 @@
enddo ! iphase
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ NGLOB_XY = NGLOB_CRUST_MANTLE
+ else
+ NGLOB_XY = 1
+ endif
+
! updates (only) acceleration w/ rotation in the crust/mantle region (touches oceans)
if(.NOT. GPU_MODE) then
- ! on CPU
- call compute_forces_el_update_accel(NGLOB_CRUST_MANTLE,veloc_crust_mantle,accel_crust_mantle, &
- two_omega_earth,rmass_crust_mantle)
- ! adjoint / kernel runs
- if (SIMULATION_TYPE == 3) &
- call compute_forces_el_update_accel(NGLOB_CRUST_MANTLE_ADJOINT,b_veloc_crust_mantle,b_accel_crust_mantle, &
- b_two_omega_earth,rmass_crust_mantle)
+ ! on CPU
+ call compute_forces_el_update_accel(NGLOB_CRUST_MANTLE,NGLOB_XY,veloc_crust_mantle,accel_crust_mantle,two_omega_earth, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NCHUNKS_VAL,ABSORBING_CONDITIONS)
+
+ ! adjoint / kernel runs
+ if (SIMULATION_TYPE == 3) &
+ call compute_forces_el_update_accel(NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_XY,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ b_two_omega_earth,rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NCHUNKS_VAL,ABSORBING_CONDITIONS)
+
else
- ! on GPU
- call kernel_3_a_cuda(Mesh_pointer, &
- deltatover2,SIMULATION_TYPE,b_deltatover2,OCEANS_VAL)
+ ! on GPU
+ call kernel_3_a_cuda(Mesh_pointer, &
+ deltatover2,SIMULATION_TYPE,b_deltatover2,OCEANS_VAL,NCHUNKS_VAL)
endif
! couples ocean with crust mantle
@@ -512,14 +522,16 @@
if(.NOT. GPU_MODE) then
! on CPU
call compute_coupling_ocean(accel_crust_mantle,b_accel_crust_mantle, &
- rmass_crust_mantle,rmass_ocean_load,normal_top_crust_mantle, &
+ rmassx_crust_mantle, rmassy_crust_mantle, rmassz_crust_mantle, &
+ rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load, &
- SIMULATION_TYPE,NSPEC2D_TOP(IREGION_CRUST_MANTLE))
+ updated_dof_ocean_load,NGLOB_XY, &
+ SIMULATION_TYPE,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+ ABSORBING_CONDITIONS)
else
! on GPU
- call compute_coupling_ocean_cuda(Mesh_pointer)
+ call compute_coupling_ocean_cuda(Mesh_pointer,NCHUNKS_VAL)
endif
endif
@@ -547,8 +559,9 @@
!=====================================================================
- subroutine compute_forces_el_update_accel(NGLOB,veloc_crust_mantle,accel_crust_mantle, &
- two_omega_earth,rmass_crust_mantle)
+ subroutine compute_forces_el_update_accel(NGLOB,NGLOB_XY,veloc_crust_mantle,accel_crust_mantle,two_omega_earth, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NCHUNKS_VAL,ABSORBING_CONDITIONS)
use constants_solver,only: CUSTOM_REAL,NDIM
@@ -558,69 +571,133 @@
implicit none
- integer :: NGLOB
+ integer :: NGLOB,NGLOB_XY,NCHUNKS_VAL
! velocity & acceleration
! crust/mantle region
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc_crust_mantle,accel_crust_mantle
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmass_crust_mantle
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmassz_crust_mantle
real(kind=CUSTOM_REAL) :: two_omega_earth
+ logical :: ABSORBING_CONDITIONS
+
! local parameters
integer :: i
! updates acceleration w/ rotation in crust/mantle region only
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+
#ifdef _HANDOPT_NEWMARK
-! way 2:
- if(imodulo_NGLOB_CRUST_MANTLE4 >= 1) then
- do i=1,imodulo_NGLOB_CRUST_MANTLE4
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ ! way 2:
+ if(imodulo_NGLOB_CRUST_MANTLE4 >= 1) then
+ do i=1,imodulo_NGLOB_CRUST_MANTLE4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
+ do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB,4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
- enddo
- endif
- do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB,4
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+
+ accel_crust_mantle(1,i+1) = accel_crust_mantle(1,i+1)*rmassx_crust_mantle(i+1) &
+ + two_omega_earth*veloc_crust_mantle(2,i+1)
+ accel_crust_mantle(2,i+1) = accel_crust_mantle(2,i+1)*rmassy_crust_mantle(i+1) &
+ - two_omega_earth*veloc_crust_mantle(1,i+1)
+ accel_crust_mantle(3,i+1) = accel_crust_mantle(3,i+1)*rmassz_crust_mantle(i+1)
+
+ accel_crust_mantle(1,i+2) = accel_crust_mantle(1,i+2)*rmassx_crust_mantle(i+2) &
+ + two_omega_earth*veloc_crust_mantle(2,i+2)
+ accel_crust_mantle(2,i+2) = accel_crust_mantle(2,i+2)*rmassy_crust_mantle(i+2) &
+ - two_omega_earth*veloc_crust_mantle(1,i+2)
+ accel_crust_mantle(3,i+2) = accel_crust_mantle(3,i+2)*rmassz_crust_mantle(i+2)
+
+ accel_crust_mantle(1,i+3) = accel_crust_mantle(1,i+3)*rmassx_crust_mantle(i+3) &
+ + two_omega_earth*veloc_crust_mantle(2,i+3)
+ accel_crust_mantle(2,i+3) = accel_crust_mantle(2,i+3)*rmassy_crust_mantle(i+3) &
+ - two_omega_earth*veloc_crust_mantle(1,i+3)
+ accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmassz_crust_mantle(i+3)
+ enddo
+#else
+ ! way 1:
+ do i=1,NGLOB
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+#endif
+
+ else
- accel_crust_mantle(1,i+1) = accel_crust_mantle(1,i+1)*rmass_crust_mantle(i+1) &
+#ifdef _HANDOPT_NEWMARK
+ ! way 2:
+ if(imodulo_NGLOB_CRUST_MANTLE4 >= 1) then
+ do i=1,imodulo_NGLOB_CRUST_MANTLE4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
+ do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB,4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+
+ accel_crust_mantle(1,i+1) = accel_crust_mantle(1,i+1)*rmassz_crust_mantle(i+1) &
+ two_omega_earth*veloc_crust_mantle(2,i+1)
- accel_crust_mantle(2,i+1) = accel_crust_mantle(2,i+1)*rmass_crust_mantle(i+1) &
+ accel_crust_mantle(2,i+1) = accel_crust_mantle(2,i+1)*rmassz_crust_mantle(i+1) &
- two_omega_earth*veloc_crust_mantle(1,i+1)
- accel_crust_mantle(3,i+1) = accel_crust_mantle(3,i+1)*rmass_crust_mantle(i+1)
+ accel_crust_mantle(3,i+1) = accel_crust_mantle(3,i+1)*rmassz_crust_mantle(i+1)
- accel_crust_mantle(1,i+2) = accel_crust_mantle(1,i+2)*rmass_crust_mantle(i+2) &
+ accel_crust_mantle(1,i+2) = accel_crust_mantle(1,i+2)*rmassz_crust_mantle(i+2) &
+ two_omega_earth*veloc_crust_mantle(2,i+2)
- accel_crust_mantle(2,i+2) = accel_crust_mantle(2,i+2)*rmass_crust_mantle(i+2) &
+ accel_crust_mantle(2,i+2) = accel_crust_mantle(2,i+2)*rmassz_crust_mantle(i+2) &
- two_omega_earth*veloc_crust_mantle(1,i+2)
- accel_crust_mantle(3,i+2) = accel_crust_mantle(3,i+2)*rmass_crust_mantle(i+2)
+ accel_crust_mantle(3,i+2) = accel_crust_mantle(3,i+2)*rmassz_crust_mantle(i+2)
- accel_crust_mantle(1,i+3) = accel_crust_mantle(1,i+3)*rmass_crust_mantle(i+3) &
+ accel_crust_mantle(1,i+3) = accel_crust_mantle(1,i+3)*rmassz_crust_mantle(i+3) &
+ two_omega_earth*veloc_crust_mantle(2,i+3)
- accel_crust_mantle(2,i+3) = accel_crust_mantle(2,i+3)*rmass_crust_mantle(i+3) &
+ accel_crust_mantle(2,i+3) = accel_crust_mantle(2,i+3)*rmassz_crust_mantle(i+3) &
- two_omega_earth*veloc_crust_mantle(1,i+3)
- accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmass_crust_mantle(i+3)
- enddo
+ accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmassz_crust_mantle(i+3)
+ enddo
#else
-! way 1:
- do i=1,NGLOB
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ ! way 1:
+ do i=1,NGLOB
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
- two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
- enddo
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
#endif
+ endif
+
end subroutine compute_forces_el_update_accel
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -212,7 +212,7 @@
tempz3l = tempz3l + displ_inner_core(3,iglob)*hp3
enddo
- if(ATTENUATION_VAL) then
+ if( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
! temporary variables used for fixing attenuation in a consistent way
tempx1l_att = 0._CUSTOM_REAL
@@ -317,7 +317,7 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- if(ATTENUATION_VAL) then
+ if( ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN ) then
! 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
@@ -337,19 +337,17 @@
duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
- epsilondev_loc(1,i,j,k) = duxdxl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(2,i,j,k) = duydyl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
endif
+ epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilondev_loc(1,i,j,k) = duxdxl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(2,i,j,k) = duydyl_att - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
else
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -276,7 +276,7 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
if(ATTENUATION_NEW_VAL) then
! takes new routines
! use first order Taylor expansion of displacement for local storage of stresses
@@ -290,72 +290,72 @@
! since we know that NGLLX = 5, this should help pipelining
iglobv5(:) = ibool(:,j,k,ispec)
- dummyx_loc_att(1,j,k) = displ_inner_core(1,iglobv5(1)) + deltat*veloc_inner_core(1,iglobv5(1))
- dummyy_loc_att(1,j,k) = displ_inner_core(2,iglobv5(1)) + deltat*veloc_inner_core(2,iglobv5(1))
- dummyz_loc_att(1,j,k) = displ_inner_core(3,iglobv5(1)) + deltat*veloc_inner_core(3,iglobv5(1))
+ dummyx_loc_att(1,j,k) = deltat*veloc_inner_core(1,iglobv5(1))
+ dummyy_loc_att(1,j,k) = deltat*veloc_inner_core(2,iglobv5(1))
+ dummyz_loc_att(1,j,k) = deltat*veloc_inner_core(3,iglobv5(1))
- dummyx_loc_att(2,j,k) = displ_inner_core(1,iglobv5(2)) + deltat*veloc_inner_core(1,iglobv5(2))
- dummyy_loc_att(2,j,k) = displ_inner_core(2,iglobv5(2)) + deltat*veloc_inner_core(2,iglobv5(2))
- dummyz_loc_att(2,j,k) = displ_inner_core(3,iglobv5(2)) + deltat*veloc_inner_core(3,iglobv5(2))
+ dummyx_loc_att(2,j,k) = deltat*veloc_inner_core(1,iglobv5(2))
+ dummyy_loc_att(2,j,k) = deltat*veloc_inner_core(2,iglobv5(2))
+ dummyz_loc_att(2,j,k) = deltat*veloc_inner_core(3,iglobv5(2))
- dummyx_loc_att(3,j,k) = displ_inner_core(1,iglobv5(3)) + deltat*veloc_inner_core(1,iglobv5(3))
- dummyy_loc_att(3,j,k) = displ_inner_core(2,iglobv5(3)) + deltat*veloc_inner_core(2,iglobv5(3))
- dummyz_loc_att(3,j,k) = displ_inner_core(3,iglobv5(3)) + deltat*veloc_inner_core(3,iglobv5(3))
+ dummyx_loc_att(3,j,k) = deltat*veloc_inner_core(1,iglobv5(3))
+ dummyy_loc_att(3,j,k) = deltat*veloc_inner_core(2,iglobv5(3))
+ dummyz_loc_att(3,j,k) = deltat*veloc_inner_core(3,iglobv5(3))
- dummyx_loc_att(4,j,k) = displ_inner_core(1,iglobv5(4)) + deltat*veloc_inner_core(1,iglobv5(4))
- dummyy_loc_att(4,j,k) = displ_inner_core(2,iglobv5(4)) + deltat*veloc_inner_core(2,iglobv5(4))
- dummyz_loc_att(4,j,k) = displ_inner_core(3,iglobv5(4)) + deltat*veloc_inner_core(3,iglobv5(4))
+ dummyx_loc_att(4,j,k) = deltat*veloc_inner_core(1,iglobv5(4))
+ dummyy_loc_att(4,j,k) = deltat*veloc_inner_core(2,iglobv5(4))
+ dummyz_loc_att(4,j,k) = deltat*veloc_inner_core(3,iglobv5(4))
- dummyx_loc_att(5,j,k) = displ_inner_core(1,iglobv5(5)) + deltat*veloc_inner_core(1,iglobv5(5))
- dummyy_loc_att(5,j,k) = displ_inner_core(2,iglobv5(5)) + deltat*veloc_inner_core(2,iglobv5(5))
- dummyz_loc_att(5,j,k) = displ_inner_core(3,iglobv5(5)) + deltat*veloc_inner_core(3,iglobv5(5))
+ dummyx_loc_att(5,j,k) = deltat*veloc_inner_core(1,iglobv5(5))
+ dummyy_loc_att(5,j,k) = deltat*veloc_inner_core(2,iglobv5(5))
+ dummyz_loc_att(5,j,k) = deltat*veloc_inner_core(3,iglobv5(5))
#else
! way 1:
do i=1,NGLLX
iglob1 = ibool(i,j,k,ispec)
- dummyx_loc_att(i,j,k) = displ_inner_core(1,iglob1) + deltat*veloc_inner_core(1,iglob1)
- dummyy_loc_att(i,j,k) = displ_inner_core(2,iglob1) + deltat*veloc_inner_core(2,iglob1)
- dummyz_loc_att(i,j,k) = displ_inner_core(3,iglob1) + deltat*veloc_inner_core(3,iglob1)
+ dummyx_loc_att(i,j,k) = deltat*veloc_inner_core(1,iglob1)
+ dummyy_loc_att(i,j,k) = deltat*veloc_inner_core(2,iglob1)
+ dummyz_loc_att(i,j,k) = deltat*veloc_inner_core(3,iglob1)
enddo
#endif
enddo
enddo
- endif
- else
- ! takes old routines
- do k=1,NGLLZ
- do j=1,NGLLY
+ else
+ ! takes old routines
+ do k=1,NGLLZ
+ do j=1,NGLLY
#ifdef _HANDOPT
- dummyx_loc_att(1,j,k) = dummyx_loc(1,j,k)
- dummyy_loc_att(1,j,k) = dummyx_loc(1,j,k)
- dummyz_loc_att(1,j,k) = dummyx_loc(1,j,k)
+ dummyx_loc_att(1,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(1,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(1,j,k) = 0._CUSTOM_REAL
- dummyx_loc_att(2,j,k) = dummyx_loc(2,j,k)
- dummyy_loc_att(2,j,k) = dummyx_loc(2,j,k)
- dummyz_loc_att(2,j,k) = dummyx_loc(2,j,k)
+ dummyx_loc_att(2,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(2,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(2,j,k) = 0._CUSTOM_REAL
- dummyx_loc_att(3,j,k) = dummyx_loc(3,j,k)
- dummyy_loc_att(3,j,k) = dummyx_loc(3,j,k)
- dummyz_loc_att(3,j,k) = dummyx_loc(3,j,k)
+ dummyx_loc_att(3,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(3,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(3,j,k) = 0._CUSTOM_REAL
- dummyx_loc_att(4,j,k) = dummyx_loc(4,j,k)
- dummyy_loc_att(4,j,k) = dummyx_loc(4,j,k)
- dummyz_loc_att(4,j,k) = dummyx_loc(4,j,k)
+ dummyx_loc_att(4,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(4,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(4,j,k) = 0._CUSTOM_REAL
- dummyx_loc_att(5,j,k) = dummyx_loc(5,j,k)
- dummyy_loc_att(5,j,k) = dummyx_loc(5,j,k)
- dummyz_loc_att(5,j,k) = dummyx_loc(5,j,k)
+ dummyx_loc_att(5,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(5,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(5,j,k) = 0._CUSTOM_REAL
#else
- do i=1,NGLLX
- dummyx_loc_att(i,j,k) = dummyx_loc(i,j,k)
- dummyy_loc_att(i,j,k) = dummyy_loc(i,j,k)
- dummyz_loc_att(i,j,k) = dummyz_loc(i,j,k)
+ do i=1,NGLLX
+ dummyx_loc_att(i,j,k) = 0._CUSTOM_REAL
+ dummyy_loc_att(i,j,k) = 0._CUSTOM_REAL
+ dummyz_loc_att(i,j,k) = 0._CUSTOM_REAL
+ enddo
+#endif
enddo
-#endif
enddo
- enddo
+ endif
endif
do j=1,m2
@@ -380,29 +380,36 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
do j=1,m2
do i=1,m1
- C1_m1_m2_5points_att(i,j) = hprime_xx(i,1)*B1_m1_m2_5points_att(1,j) + &
+ C1_m1_m2_5points_att(i,j) = C1_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B1_m1_m2_5points_att(1,j) + &
hprime_xx(i,2)*B1_m1_m2_5points_att(2,j) + &
hprime_xx(i,3)*B1_m1_m2_5points_att(3,j) + &
hprime_xx(i,4)*B1_m1_m2_5points_att(4,j) + &
hprime_xx(i,5)*B1_m1_m2_5points_att(5,j)
- C2_m1_m2_5points_att(i,j) = hprime_xx(i,1)*B2_m1_m2_5points_att(1,j) + &
+ C2_m1_m2_5points_att(i,j) = C2_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B2_m1_m2_5points_att(1,j) + &
hprime_xx(i,2)*B2_m1_m2_5points_att(2,j) + &
hprime_xx(i,3)*B2_m1_m2_5points_att(3,j) + &
hprime_xx(i,4)*B2_m1_m2_5points_att(4,j) + &
hprime_xx(i,5)*B2_m1_m2_5points_att(5,j)
- C3_m1_m2_5points_att(i,j) = hprime_xx(i,1)*B3_m1_m2_5points_att(1,j) + &
+ C3_m1_m2_5points_att(i,j) = C3_m1_m2_5points(i,j) + &
+ hprime_xx(i,1)*B3_m1_m2_5points_att(1,j) + &
hprime_xx(i,2)*B3_m1_m2_5points_att(2,j) + &
hprime_xx(i,3)*B3_m1_m2_5points_att(3,j) + &
hprime_xx(i,4)*B3_m1_m2_5points_att(4,j) + &
hprime_xx(i,5)*B3_m1_m2_5points_att(5,j)
enddo
enddo
+ else
+ tempx1_att(:,:,:) = 0._CUSTOM_REAL
+ tempy1_att(:,:,:) = 0._CUSTOM_REAL
+ tempz1_att(:,:,:) = 0._CUSTOM_REAL
endif
do j=1,m1
@@ -430,25 +437,28 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
do j=1,m1
do i=1,m1
! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
do k = 1,NGLLX
- tempx2_att(i,j,k) = dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ tempx2_att(i,j,k) = tempx2(i,j,k) + &
+ dummyx_loc_att(i,1,k)*hprime_xxT(1,j) + &
dummyx_loc_att(i,2,k)*hprime_xxT(2,j) + &
dummyx_loc_att(i,3,k)*hprime_xxT(3,j) + &
dummyx_loc_att(i,4,k)*hprime_xxT(4,j) + &
dummyx_loc_att(i,5,k)*hprime_xxT(5,j)
- tempy2_att(i,j,k) = dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ tempy2_att(i,j,k) = tempy2(i,j,k) + &
+ dummyy_loc_att(i,1,k)*hprime_xxT(1,j) + &
dummyy_loc_att(i,2,k)*hprime_xxT(2,j) + &
dummyy_loc_att(i,3,k)*hprime_xxT(3,j) + &
dummyy_loc_att(i,4,k)*hprime_xxT(4,j) + &
dummyy_loc_att(i,5,k)*hprime_xxT(5,j)
- tempz2_att(i,j,k) = dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
+ tempz2_att(i,j,k) = tempz2(i,j,k) + &
+ dummyz_loc_att(i,1,k)*hprime_xxT(1,j) + &
dummyz_loc_att(i,2,k)*hprime_xxT(2,j) + &
dummyz_loc_att(i,3,k)*hprime_xxT(3,j) + &
dummyz_loc_att(i,4,k)*hprime_xxT(4,j) + &
@@ -456,6 +466,10 @@
enddo
enddo
enddo
+ else
+ tempx2_att(:,:,:) = 0._CUSTOM_REAL
+ tempy2_att(:,:,:) = 0._CUSTOM_REAL
+ tempz2_att(:,:,:) = 0._CUSTOM_REAL
endif
do j=1,m1
@@ -480,29 +494,36 @@
enddo
enddo
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
do j=1,m1
do i=1,m2
- C1_mxm_m2_m1_5points_att(i,j) = A1_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ C1_mxm_m2_m1_5points_att(i,j) = C1_mxm_m2_m1_5points(i,j) + &
+ A1_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
A1_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
A1_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
A1_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
A1_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
- C2_mxm_m2_m1_5points_att(i,j) = A2_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ C2_mxm_m2_m1_5points_att(i,j) = C2_mxm_m2_m1_5points(i,j) + &
+ A2_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
A2_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
A2_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
A2_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
A2_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
- C3_mxm_m2_m1_5points_att(i,j) = A3_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
+ C3_mxm_m2_m1_5points_att(i,j) = C3_mxm_m2_m1_5points(i,j) + &
+ A3_mxm_m2_m1_5points_att(i,1)*hprime_xxT(1,j) + &
A3_mxm_m2_m1_5points_att(i,2)*hprime_xxT(2,j) + &
A3_mxm_m2_m1_5points_att(i,3)*hprime_xxT(3,j) + &
A3_mxm_m2_m1_5points_att(i,4)*hprime_xxT(4,j) + &
A3_mxm_m2_m1_5points_att(i,5)*hprime_xxT(5,j)
enddo
enddo
+ else
+ tempx3_att(:,:,:) = 0._CUSTOM_REAL
+ tempy3_att(:,:,:) = 0._CUSTOM_REAL
+ tempz3_att(:,:,:) = 0._CUSTOM_REAL
endif
do k=1,NGLLZ
@@ -545,7 +566,7 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- if(ATTENUATION_VAL) then
+ if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
@@ -565,20 +586,18 @@
duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att
! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- epsilondev_loc(1,i,j,k) = duxdxl_att - templ
- epsilondev_loc(2,i,j,k) = duydyl_att - templ
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+ if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
endif
+ templ = ONE_THIRD * (duxdxl_att + duydyl_att + duzdzl_att)
+ epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilondev_loc(1,i,j,k) = duxdxl_att - templ
+ epsilondev_loc(2,i,j,k) = duydyl_att - templ
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl_att
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
else
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_crust_mantle.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_crust_mantle.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -84,7 +84,7 @@
call read_abs(0,absorb_xmin_crust_mantle,reclen_xmin_crust_mantle,NSTEP-it+1)
endif
- if( .NOT. GPU_MODE) then
+ if ( .NOT. GPU_MODE) then
! on CPU
do ispec2D=1,nspec2D_xmin_crust_mantle
@@ -158,7 +158,7 @@
ispec=ibelm_xmax_crust_mantle(ispec2D)
- ! exclude elements that are not on absorbing edges
+ ! exclude elements that are not on absorbing edges
if(nkmin_xi_crust_mantle(2,ispec2D) == 0 .or. njmin_crust_mantle(2,ispec2D) == 0) cycle
i=NGLLX
@@ -294,7 +294,7 @@
ispec=ibelm_ymax_crust_mantle(ispec2D)
- ! exclude elements that are not on absorbing edges
+ ! exclude elements that are not on absorbing edges
if(nkmin_eta_crust_mantle(2,ispec2D) == 0 .or. nimin_crust_mantle(2,ispec2D) == 0) cycle
j=NGLLY
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -140,6 +140,16 @@
! frees dynamically allocated memory
+ ! mass matrices
+
+ deallocate(rmassx_crust_mantle)
+ deallocate(rmassy_crust_mantle)
+ deallocate(rmassz_crust_mantle)
+
+ deallocate(rmass_outer_core)
+ deallocate(rmass_inner_core)
+
+
! mpi buffers
deallocate(buffer_send_vector_crust_mantle,buffer_recv_vector_crust_mantle, &
request_send_vector_crust_mantle,request_recv_vector_crust_mantle)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -573,6 +573,6 @@
endif
! frees allocated memory on GPU
- call prepare_cleanup_device(Mesh_pointer)
+ call prepare_cleanup_device(Mesh_pointer,NCHUNKS_VAL)
end subroutine it_transfer_from_GPU
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-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -193,17 +193,32 @@
if(minval(rmass_ocean_load) <= 0._CUSTOM_REAL) &
call exit_MPI(myrank,'negative mass matrix term for the oceans')
endif
- if(minval(rmass_crust_mantle) <= 0._CUSTOM_REAL) &
- call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
+
+ ! add C*deltat/2 contribution to the mass matrices on Stacey edges
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if(minval(rmassx_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
+ if(minval(rmassy_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
+ endif
+
+ if(minval(rmassz_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
if(minval(rmass_inner_core) <= 0._CUSTOM_REAL) &
- call exit_MPI(myrank,'negative mass matrix term for the inner core')
+ call exit_MPI(myrank,'negative mass matrix term for the inner core')
if(minval(rmass_outer_core) <= 0._CUSTOM_REAL) &
- call exit_MPI(myrank,'negative mass matrix term for the outer core')
+ call exit_MPI(myrank,'negative mass matrix term for the outer core')
! for efficiency, invert final mass matrix once and for all on each slice
if(OCEANS_VAL) rmass_ocean_load = 1._CUSTOM_REAL / rmass_ocean_load
- rmass_crust_mantle = 1._CUSTOM_REAL / rmass_crust_mantle
+ ! add C*deltat/2 contribution to the mass matrices on Stacey edges
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ rmassx_crust_mantle = 1._CUSTOM_REAL / rmassx_crust_mantle
+ rmassy_crust_mantle = 1._CUSTOM_REAL / rmassy_crust_mantle
+ endif
+
+ rmassz_crust_mantle = 1._CUSTOM_REAL / rmassz_crust_mantle
rmass_outer_core = 1._CUSTOM_REAL / rmass_outer_core
rmass_inner_core = 1._CUSTOM_REAL / rmass_inner_core
@@ -233,8 +248,22 @@
endif
! crust and mantle
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ call assemble_MPI_scalar_ext_mesh(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, &
+ rmassx_crust_mantle, &
+ num_interfaces_crust_mantle,max_nibool_interfaces_crust_mantle, &
+ nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle,&
+ my_neighbours_crust_mantle)
+
+ call assemble_MPI_scalar_ext_mesh(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, &
+ rmassy_crust_mantle, &
+ num_interfaces_crust_mantle,max_nibool_interfaces_crust_mantle, &
+ nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle,&
+ my_neighbours_crust_mantle)
+ endif
+
call assemble_MPI_scalar_ext_mesh(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, &
- rmass_crust_mantle, &
+ rmassz_crust_mantle, &
num_interfaces_crust_mantle,max_nibool_interfaces_crust_mantle, &
nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle,&
my_neighbours_crust_mantle)
@@ -1265,32 +1294,36 @@
! crust/mantle region
if(myrank == 0 ) write(IMAIN,*) " loading crust/mantle region"
+
call prepare_crust_mantle_device(Mesh_pointer, &
- xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
- etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
- rhostore_crust_mantle, &
- kappavstore_crust_mantle,muvstore_crust_mantle, &
- kappahstore_crust_mantle,muhstore_crust_mantle, &
- eta_anisostore_crust_mantle, &
- rmass_crust_mantle, &
- normal_top_crust_mantle, &
- ibelm_top_crust_mantle, &
- ibelm_bottom_crust_mantle, &
- ibool_crust_mantle, &
- xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
- ispec_is_tiso_crust_mantle, &
- c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
- c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle, &
- c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle, &
- c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle, &
- c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
- c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
- c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
- num_phase_ispec_crust_mantle,phase_ispec_inner_crust_mantle, &
- nspec_outer_crust_mantle,nspec_inner_crust_mantle, &
- NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE))
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
+ rhostore_crust_mantle, &
+ kappavstore_crust_mantle,muvstore_crust_mantle, &
+ kappahstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle, &
+ rmassx_crust_mantle, &
+ rmassy_crust_mantle, &
+ rmassz_crust_mantle, &
+ normal_top_crust_mantle, &
+ ibelm_top_crust_mantle, &
+ ibelm_bottom_crust_mantle, &
+ ibool_crust_mantle, &
+ xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+ ispec_is_tiso_crust_mantle, &
+ c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
+ c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle, &
+ c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle, &
+ c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle, &
+ c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
+ c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
+ c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
+ num_phase_ispec_crust_mantle,phase_ispec_inner_crust_mantle, &
+ nspec_outer_crust_mantle,nspec_inner_crust_mantle, &
+ NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+ NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE), &
+ NCHUNKS_VAL)
call sync_all()
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -35,8 +35,8 @@
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,ispec_is_tiso, &
- rmass,rmass_ocean_load,nspec,nglob, &
+ ibool,idoubling,ispec_is_tiso,nglob_xy,nglob, &
+ rmassx,rmassy,rmassz,rmass_ocean_load,nspec, &
READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY, &
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS,LOCAL_PATH,ABSORBING_CONDITIONS)
@@ -54,7 +54,7 @@
character(len=150) LOCAL_PATH
- integer nspec,nglob
+ integer nspec,nglob,nglob_xy
integer nspec_iso,nspec_tiso,nspec_ani
@@ -83,9 +83,12 @@
real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,nspec)
real(kind=CUSTOM_REAL) rho_vs(NGLLX,NGLLY,NGLLZ,nspec)
- ! mass matrix and additional ocean load mass matrix
- real(kind=CUSTOM_REAL), dimension(nglob) :: rmass,rmass_ocean_load
+ ! mass matrices and additional ocean load mass matrix
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
+
! global addressing
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
@@ -173,9 +176,21 @@
endif
- ! mass matrix
- read(IIN) rmass
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) then
+ read(IIN) rmassx
+ read(IIN) rmassy
+ endif
+ read(IIN) rmassz
+
! read additional ocean load mass matrix
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) read(IIN) rmass_ocean_load
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -90,7 +90,7 @@
implicit none
! local parameters
- integer :: nspec_iso,nspec_tiso,nspec_ani
+ integer :: nspec_iso,nspec_tiso,nspec_ani,NGLOB_XY
logical :: READ_KAPPA_MU,READ_TISO
! dummy array that does not need to be actually read
integer, dimension(:),allocatable :: dummy_i
@@ -122,6 +122,27 @@
! sets number of top elements for surface movies & noise tomography
NSPEC_TOP = NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ ! allocates mass matrices in this slice (will be fully assembled in the solver)
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ NGLOB_XY = NGLOB_CRUST_MANTLE
+ else
+ NGLOB_XY = 1
+ endif
+
+ allocate(rmassx_crust_mantle(NGLOB_XY),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+ allocate(rmassy_crust_mantle(NGLOB_XY),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+ allocate(rmassz_crust_mantle(NGLOB_CRUST_MANTLE),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+
! reads databases file
call read_arrays_solver(IREGION_CRUST_MANTLE,myrank, &
rho_vp_crust_mantle,rho_vs_crust_mantle, &
@@ -139,10 +160,8 @@
c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
- ibool_crust_mantle,dummy_i, &
- ispec_is_tiso_crust_mantle, &
- rmass_crust_mantle,rmass_ocean_load, &
- NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE, &
+ ibool_crust_mantle,dummy_i,ispec_is_tiso_crust_mantle,NGLOB_XY,NGLOB_CRUST_MANTLE, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle,rmass_ocean_load,NSPEC_CRUST_MANTLE, &
READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS)
@@ -168,14 +187,16 @@
implicit none
! local parameters
- integer :: nspec_iso,nspec_tiso,nspec_ani
+ integer :: nspec_iso,nspec_tiso,nspec_ani,NGLOB_XY
logical :: READ_KAPPA_MU,READ_TISO
+ integer :: ier
! dummy array that does not need to be actually read
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,1) :: dummy_array
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: dummy_rmass
- logical, dimension(:),allocatable:: dummy_ispec_is_tiso
- integer, dimension(:),allocatable :: dummy_idoubling_outer_core
+ logical, dimension(:), allocatable :: dummy_ispec_is_tiso
+ integer, dimension(:), allocatable :: dummy_idoubling_outer_core
! outer core (no anisotropy nor S velocity)
! rmass_ocean_load is not used in this routine because it is meaningless in the outer core
@@ -186,9 +207,23 @@
nspec_ani = 1
! dummy allocation
+ NGLOB_XY = 1
+
+ allocate(dummy_rmass(NGLOB_XY))
allocate(dummy_ispec_is_tiso(NSPEC_OUTER_CORE))
allocate(dummy_idoubling_outer_core(NSPEC_OUTER_CORE))
+ ! allocates mass matrices in this slice (will be fully assembled in the solver)
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ allocate(rmass_outer_core(NGLOB_OUTER_CORE),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+
call read_arrays_solver(IREGION_OUTER_CORE,myrank, &
vp_outer_core,dummy_array, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
@@ -205,14 +240,14 @@
dummy_array,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
- ibool_outer_core,dummy_idoubling_outer_core,dummy_ispec_is_tiso, &
- rmass_outer_core,rmass_ocean_load, &
- NSPEC_OUTER_CORE,NGLOB_OUTER_CORE, &
+ ibool_outer_core,dummy_idoubling_outer_core,dummy_ispec_is_tiso,NGLOB_XY,NGLOB_OUTER_CORE, &
+ dummy_rmass,dummy_rmass,rmass_outer_core,rmass_ocean_load,NSPEC_OUTER_CORE, &
READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS)
deallocate(dummy_idoubling_outer_core)
deallocate(dummy_ispec_is_tiso)
+ deallocate(dummy_rmass)
! check that the number of points in this slice is correct
if(minval(ibool_outer_core(:,:,:,:)) /= 1 .or. &
@@ -234,11 +269,13 @@
implicit none
! local parameters
- integer :: nspec_iso,nspec_tiso,nspec_ani
+ integer :: nspec_iso,nspec_tiso,nspec_ani,NGLOB_XY
logical :: READ_KAPPA_MU,READ_TISO
+ integer :: ier
! dummy array that does not need to be actually read
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,1) :: dummy_array
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: dummy_rmass
logical, dimension(:),allocatable:: dummy_ispec_is_tiso
! inner core (no anisotropy)
@@ -254,8 +291,22 @@
endif
! dummy allocation
+ NGLOB_XY = 1
+
+ allocate(dummy_rmass(NGLOB_XY))
allocate(dummy_ispec_is_tiso(NSPEC_INNER_CORE))
+ ! allocates mass matrices in this slice (will be fully assembled in the solver)
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ allocate(rmass_inner_core(NGLOB_INNER_CORE),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+
call read_arrays_solver(IREGION_INNER_CORE,myrank, &
dummy_array,dummy_array, &
xstore_inner_core,ystore_inner_core,zstore_inner_core, &
@@ -272,13 +323,13 @@
dummy_array,dummy_array,dummy_array, &
c44store_inner_core,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
- ibool_inner_core,idoubling_inner_core,dummy_ispec_is_tiso, &
- rmass_inner_core,rmass_ocean_load, &
- NSPEC_INNER_CORE,NGLOB_INNER_CORE, &
+ ibool_inner_core,idoubling_inner_core,dummy_ispec_is_tiso,NGLOB_XY,NGLOB_INNER_CORE, &
+ dummy_rmass,dummy_rmass,rmass_inner_core,rmass_ocean_load,NSPEC_INNER_CORE, &
READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS)
deallocate(dummy_ispec_is_tiso)
+ deallocate(dummy_rmass)
! check that the number of points in this slice is correct
if(minval(ibool_inner_core(:,:,:,:)) /= 1 .or. maxval(ibool_inner_core(:,:,:,:)) /= NGLOB_INNER_CORE) &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2012-06-15 01:06:32 UTC (rev 20374)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2012-06-15 11:07:03 UTC (rev 20375)
@@ -363,8 +363,17 @@
! flag for transversely isotropic elements
logical, dimension(NSPEC_CRUST_MANTLE) :: ispec_is_tiso_crust_mantle
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassz_crust_mantle
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
@@ -509,7 +518,7 @@
integer, dimension(NSPEC_INNER_CORE) :: idoubling_inner_core
! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: rmass_inner_core
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inner_core
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: &
@@ -597,7 +606,7 @@
rhostore_outer_core,kappavstore_outer_core
! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: rmass_outer_core
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_outer_core
! velocity potential
real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: &
More information about the CIG-COMMITS
mailing list