[cig-commits] r19214 - in seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src: cuda shared specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Nov 17 19:01:34 PST 2011
Author: danielpeter
Date: 2011-11-17 19:01:34 -0800 (Thu, 17 Nov 2011)
New Revision: 19214
Modified:
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_elastic_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
Log:
updates arrays for gravity
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_elastic_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_elastic_cuda.cu 2011-11-17 23:51:03 UTC (rev 19213)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_elastic_cuda.cu 2011-11-18 03:01:34 UTC (rev 19214)
@@ -434,8 +434,6 @@
minus_g = d_minus_g[iglob];
minus_dg = d_minus_deriv_gravity[iglob];
- rhol = d_rhostore[working_element*NGLL3_PADDED + tx];
-
// Cartesian components of the gravitational acceleration
//gxl = 0.f;
//gyl = 0.f;
@@ -451,6 +449,8 @@
//Hxzl = 0.f;
//Hyzl = 0.f;
+ rhol = d_rhostore[working_element*NGLL3_PADDED + tx];
+
// get displacement and multiply by density to compute G tensor
// G = rho [ sg - (s * g) I ]
sx_l = rhol * s_dummyx_loc[tx]; // d_displ[iglob*3];
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2011-11-17 23:51:03 UTC (rev 19213)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2011-11-18 03:01:34 UTC (rev 19214)
@@ -92,7 +92,8 @@
{
fprintf(stderr,"Error after %s: %s\n", kernel_name, cudaGetErrorString(err));
pause_for_debugger(0);
- exit(1);
+ free(kernel_name);
+ exit(EXIT_FAILURE);
}
}
@@ -105,6 +106,7 @@
#ifdef USE_MPI
MPI_Abort(MPI_COMM_WORLD,1);
#endif
+ free(info);
exit(EXIT_FAILURE);
return;
}
@@ -119,8 +121,8 @@
fflush(stdout);
#ifdef USE_MPI
MPI_Abort(MPI_COMM_WORLD,1);
-#endif
- exit(0);
+#endif
+ exit(EXIT_FAILURE);
}
return;
}
@@ -135,7 +137,7 @@
cudaError_t cuda_status = cudaMemGetInfo( &free_byte, &total_byte ) ;
if ( cudaSuccess != cuda_status ){
printf("Error: cudaMemGetInfo fails, %s \n", cudaGetErrorString(cuda_status) );
- exit(1);
+ exit(EXIT_FAILURE);
}
*free_db = (double)free_byte ;
@@ -156,7 +158,7 @@
get_free_memory(&free_db,&used_db,&total_db);
- sprintf(filename,"../in_out_files/OUTPUT_FILES/gpu_mem_usage_proc_%06d.txt",myrank);
+ sprintf(filename,"../in_out_files/OUTPUT_FILES/gpu_device_mem_usage_proc_%06d.txt",myrank);
fp = fopen(filename,"a+");
fprintf(fp,"%d: @%s GPU memory usage: used = %f MB, free = %f MB, total = %f MB\n", myrank, info_str,
used_db/1024.0/1024.0, free_db/1024.0/1024.0, total_db/1024.0/1024.0);
@@ -758,6 +760,9 @@
mp->nspec_acoustic = *nspec_acoustic;
mp->nspec_elastic = *nspec_elastic;
+ // gravity flag initialization
+ mp->gravity = 0;
+
#ifdef ENABLE_VERY_SLOW_ERROR_CHECKING
exit_on_cuda_error("prepare_constants_device");
#endif
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90 2011-11-17 23:51:03 UTC (rev 19213)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90 2011-11-18 03:01:34 UTC (rev 19214)
@@ -166,15 +166,15 @@
call get_value_string(CMTSOLUTION, 'solver.CMTSOLUTION',&
IN_DATA_FILES_PATH(1:len_trim(IN_DATA_FILES_PATH))//'CMTSOLUTION')
- open(unit=1,file=CMTSOLUTION,iostat=ios,status='old',action='read')
+ open(unit=21,file=trim(CMTSOLUTION),iostat=ios,status='old',action='read')
if(ios /= 0) stop 'error opening CMTSOLUTION file'
icounter = 0
do while(ios == 0)
- read(1,"(a)",iostat=ios) dummystring
+ read(21,"(a)",iostat=ios) dummystring
if(ios == 0) icounter = icounter + 1
enddo
- close(1)
+ close(21)
if(mod(icounter,NLINES_PER_CMTSOLUTION_SOURCE) /= 0) &
stop 'total number of lines in CMTSOLUTION file should be a multiple of NLINES_PER_CMTSOLUTION_SOURCE'
@@ -183,27 +183,27 @@
if(NSOURCES < 1) stop 'need at least one source in CMTSOLUTION file'
! compute the minimum value of hdur in CMTSOLUTION file
- open(unit=1,file=CMTSOLUTION,status='old',action='read')
+ open(unit=21,file=trim(CMTSOLUTION),status='old',action='read')
minval_hdur = HUGEVAL
do isource = 1,NSOURCES
! skip other information
do idummy = 1,3
- read(1,"(a)") dummystring
+ read(21,"(a)") dummystring
enddo
! read half duration and compute minimum
- read(1,"(a)") dummystring
+ read(21,"(a)") dummystring
read(dummystring(15:len_trim(dummystring)),*) hdur
minval_hdur = min(minval_hdur,hdur)
! skip other information
do idummy = 1,9
- read(1,"(a)") dummystring
+ read(21,"(a)") dummystring
enddo
enddo
- close(1)
+ close(21)
! one cannot use a Heaviside source for the movies
if((MOVIE_SURFACE .or. MOVIE_VOLUME) .and. sqrt(minval_hdur**2 + HDUR_MOVIE**2) < TINYVAL) &
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2011-11-17 23:51:03 UTC (rev 19213)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2011-11-18 03:01:34 UTC (rev 19214)
@@ -477,19 +477,21 @@
RICB,RCMB,RTOPDDOUBLEPRIME, &
R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
+ dg = 4.0d0*rho - 2.0d0*g/radius
+
! re-dimensionalize
- radius = radius * R_EARTH ! in m
- vp = vp * R_EARTH*dsqrt(PI*GRAV*RHOAV) ! in m / s
- rho = rho * RHOAV ! in kg / m^3
g = g * R_EARTH*(PI*GRAV*RHOAV) ! in m / s^2 ( should be around 10 m/s^2 )
+ dg = dg * R_EARTH*(PI*GRAV*RHOAV) / R_EARTH ! gradient d/dz g , in 1/s^2
- dg = 4.0d0*rho - 2.0d0*g/radius
-
minus_deriv_gravity(iglob) = - dg
- minus_g(iglob) = - g ! in negative z-direction - g ! / vp**2
+ minus_g(iglob) = - g ! in negative z-direction
! debug
!if( iglob == 1 .or. iglob == 1000 .or. iglob == 10000 ) then
+ ! ! re-dimensionalize
+ ! radius = radius * R_EARTH ! in m
+ ! vp = vp * R_EARTH*dsqrt(PI*GRAV*RHOAV) ! in m / s
+ ! rho = rho * RHOAV ! in kg / m^3
! print*,'gravity: radius=',radius,'g=',g,'depth=',radius-R_EARTH
! print*,'vp=',vp,'rho=',rho,'kappa=',(vp**2) * rho
! print*,'minus_g..=',minus_g(iglob)
@@ -981,10 +983,12 @@
endif ! NOISE_TOMOGRAPHY
! prepares gravity arrays
- call prepare_fields_gravity_device(Mesh_pointer,GRAVITY, &
+ if( GRAVITY ) then
+ call prepare_fields_gravity_device(Mesh_pointer,GRAVITY, &
minus_deriv_gravity,minus_g,wgll_cube,&
ACOUSTIC_SIMULATION,rhostore)
-
+ endif
+
! sends initial data to device
! puts acoustic initial fields onto GPU
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2011-11-17 23:51:03 UTC (rev 19213)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2011-11-18 03:01:34 UTC (rev 19214)
@@ -186,6 +186,20 @@
read(27) rho_vp
read(27) rho_vs
+ ! checks if rhostore is available for gravity
+ if( GRAVITY ) then
+ if( .not. ACOUSTIC_SIMULATION ) then
+ ! rho array needed for gravity
+ allocate(rhostore(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rhostore'
+ ! extract rho information from mu = rho * vs * vs and rho_vs = rho * vs
+ where( mustore > TINYVAL )
+ rhostore = (rho_vs*rho_vs) / mustore
+ elsewhere
+ rhostore(:,:,:,:) = 0.0_CUSTOM_REAL
+ endwhere
+ endif
+ endif
else
! no elastic attenuation & anisotropy
ATTENUATION = .false.
More information about the CIG-COMMITS
mailing list