[cig-commits] r19406 - in seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER: doc/paper_draft examples src/cuda src/generate_databases src/shared src/specfem3D utils
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Jan 20 16:35:49 PST 2012
Author: danielpeter
Date: 2012-01-20 16:35:49 -0800 (Fri, 20 Jan 2012)
New Revision: 19406
Modified:
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/doc/paper_draft/paper_sesame_2.0.tex
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/examples/README
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/check_fields_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_acoustic_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_elastic_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_kernels_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_stacey_acoustic_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/mesh_constants_cuda.h
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/noise_tomography_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/write_seismograms_cuda.cu
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/compute_arrays_source.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/write_VTK_data.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_acoustic_el.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_ac.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_elastic.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/locate_partition.f90
Log:
updates routines and Makefiles for serial configuration (--without-mpi)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/doc/paper_draft/paper_sesame_2.0.tex
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/doc/paper_draft/paper_sesame_2.0.tex 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/doc/paper_draft/paper_sesame_2.0.tex 2012-01-21 00:35:49 UTC (rev 19406)
@@ -1022,7 +1022,7 @@
\eequ
It is this last kernel expression that is actually implemented,
since the values for $ \partial_t^2 \phi $ and $ \partial_t^2 \phi^\dag$ are obtained at each time step
-in the Newark time scheme used to propagate acoustic waves.
+in the Newmark time scheme used to propagate acoustic waves.
%% ------------------------------------------------------------------------ %%
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/examples/README
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/examples/README 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/examples/README 2012-01-21 00:35:49 UTC (rev 19406)
@@ -23,8 +23,8 @@
- noise_tomography/: uses the homogeneous_halfspace model for noise kernel simulations
+- Mount_StHelens/: a model with topography data included
-
To familiarize yourself with the new in-house mesher xmeshfem3D (no CUBIT needed),
you can follow the examples in:
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/check_fields_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/check_fields_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/check_fields_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -29,7 +29,10 @@
#include <stdio.h>
#include <cuda.h>
#include <cublas.h>
+
+#ifdef WITH_MPI
#include <mpi.h>
+#endif
#include <sys/time.h>
#include <sys/resource.h>
@@ -71,7 +74,11 @@
TRACE("check_max_norm_vector");
int procid;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
+#else
+ procid = 0;
+#endif
realw maxnorm=0;
int maxloc;
for(int i=0;i<*size;i++) {
@@ -226,7 +233,11 @@
printf("rel error = %f, maxerr = %e @ %d\n",diff2/sum,maxerr,maxerrorloc);
int myrank;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+#else
+ myrank = 0;
+#endif
if(myrank == 0) {
for(int i=maxerrorloc;i>maxerrorloc-5;i--) {
printf("[%d]: %e vs. %e\n",i,vector1[i],vector2[i]);
@@ -253,7 +264,11 @@
Mesh* mp = (Mesh*)(*Mesh_pointer);
int procid;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
+#else
+ procid = 0;
+#endif
int size = *sizef;
int it = *itf;
realw* accel_cpy = (realw*)malloc(size*sizeof(realw));
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -69,14 +69,14 @@
// if(igll<NGLL2) {
// "-1" from index values to convert from Fortran-> C indexing
- ispec = coupling_ac_el_ispec[iface]-1;
+ ispec = coupling_ac_el_ispec[iface] - 1;
if(ispec_is_inner[ispec] == phase_is_inner ) {
i = coupling_ac_el_ijk[INDEX3(NDIM,NGLL2,0,igll,iface)] - 1;
j = coupling_ac_el_ijk[INDEX3(NDIM,NGLL2,1,igll,iface)] - 1;
k = coupling_ac_el_ijk[INDEX3(NDIM,NGLL2,2,igll,iface)] - 1;
- iglob = ibool[INDEX4(5,5,5,i,j,k,ispec)]- 1;
+ iglob = ibool[INDEX4(5,5,5,i,j,k,ispec)] - 1;
// elastic displacement on global point
displ_x = displ[iglob*3] ; // (1,iglob)
@@ -97,7 +97,7 @@
// continuity of pressure and normal displacement on global point
- // note: newark time scheme together with definition of scalar potential:
+ // note: Newmark time scheme together with definition of scalar potential:
// pressure = - chi_dot_dot
// requires that this coupling term uses the updated displacement at time step [t+delta_t],
// which is done at the very beginning of the time loop
@@ -212,14 +212,14 @@
// if(igll<NGLL2) {
// "-1" from index values to convert from Fortran-> C indexing
- ispec = coupling_ac_el_ispec[iface]-1;
+ ispec = coupling_ac_el_ispec[iface] - 1;
if(ispec_is_inner[ispec] == phase_is_inner ) {
i = coupling_ac_el_ijk[INDEX3(NDIM,NGLL2,0,igll,iface)] - 1;
j = coupling_ac_el_ijk[INDEX3(NDIM,NGLL2,1,igll,iface)] - 1;
k = coupling_ac_el_ijk[INDEX3(NDIM,NGLL2,2,igll,iface)] - 1;
- iglob = ibool[INDEX4(5,5,5,i,j,k,ispec)]- 1;
+ iglob = ibool[INDEX4(5,5,5,i,j,k,ispec)] - 1;
// gets associated normal on GLL point
// note: normal points away from acoustic element
@@ -238,18 +238,28 @@
// note: uses potential chi such that displacement s = grad(chi),
// pressure becomes: p = - kappa ( div( s ) ) = rho ( - dot_dot_chi + g * s )
// g only acting in negative z-direction
- // daniel: TODO - check coupling would be displ * nz correct?
+
+ // daniel: TODO - check gravity and coupling would be displ * nz correct?
pressure = rhol*( - potential_dot_dot_acoustic[iglob]
+ minus_g[iglob] * displ[iglob*3+2] );
+ //daniel: TODO - check gravity and coupling
+ //pressure = - potential_dot_dot_acoustic[iglob] ;
+ //if( iface == 128 && igll == 5 ){
+ // printf("coupling acoustic: %f %f \n",potential_dot_dot_acoustic[iglob],
+ // minus_g[iglob] * displ[iglob*3+2]);
+ //}
+
}else{
// no gravity: uses potential chi such that displacement s = 1/rho grad(chi)
+ // pressure p = - kappa ( div( s ) ) then becomes: p = - dot_dot_chi
+ // ( multiplied with factor 1/kappa due to setup of equation of motion )
pressure = - potential_dot_dot_acoustic[iglob];
}
// continuity of displacement and pressure on global point
//
- // note: newark time scheme together with definition of scalar potential:
+ // note: Newmark time scheme together with definition of scalar potential:
// pressure = - chi_dot_dot
// requires that this coupling term uses the *UPDATED* pressure (chi_dot_dot), i.e.
// pressure at time step [t + delta_t]
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_acoustic_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_acoustic_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_acoustic_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -411,9 +411,26 @@
// gravity term: 1/kappa grad(chi) * g
// assumes that g only acts in (negative) z-direction
kappa_invl = 1.f / d_kappastore[working_element*NGLL3 + tx];
-
iglob = d_ibool[working_element*NGLL3 + tx]-1;
+
+ // daniel: TODO - check gravity
+// if( kappa_invl <= 0.0f ){
+// printf("kappa error: %f %f\n",kappa_invl,d_kappastore[working_element*NGLL3 + tx]);
+// printf("kappa error: thread %d %d \n",tx,working_element);
+// asm("trap;");
+// }
+// if( iglob <= 0 ){
+// printf("iglob error: %d %d %d \n",iglob,tx,working_element);
+// asm("trap;");
+// }
+
gravity_term = minus_g[iglob] * kappa_invl * jacobianl * wgll_cube[tx] * dpotentialdzl;
+
+ // daniel: TODO - check gravity
+ //gravity_term = 0.f;
+ //if( iglob == 5 ){
+ // printf("iglob infos: %f %f %f %f %f \n",minus_g[iglob],kappa_invl,jacobianl,wgll_cube[tx],dpotentialdzl);
+ //}
}
// density (reciproc)
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 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_forces_elastic_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -479,11 +479,17 @@
//rho_s_H1 = fac1 * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl);
//rho_s_H2 = fac1 * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl);
//rho_s_H3 = fac1 * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl);
+
// only non-zero z-direction
*rho_s_H1 = factor * sx_l * Hxxl ; // 0.f;
*rho_s_H2 = factor * sy_l * Hyyl ; // 0.f;
*rho_s_H3 = factor * sz_l * Hzzl ;
+ // debug
+ //*rho_s_H1 = 0.f;
+ //*rho_s_H2 = 0.f;
+ //*rho_s_H3 = 0.f ;
+
}
/* ----------------------------------------------------------------------------------------------- */
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_kernels_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_kernels_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_kernels_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -315,7 +315,8 @@
gammazl = d_gammaz[offset];
if( gravity ){
- rho_invl = 1.0f;
+ // daniel: TODO - check gravity case here
+ rho_invl = 1.0f / rhol;
}else{
rho_invl = 1.0f / rhol;
}
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_stacey_acoustic_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_stacey_acoustic_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/compute_stacey_acoustic_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -90,7 +90,7 @@
// velocity
if( gravity ){
- // daniel: TODO - check stacey here...
+ // daniel: TODO - check gravity and stacey condition here...
// uses a potential definition of: s = grad(chi)
vel = potential_dot_acoustic[iglob] / rhol ;
}else{
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/mesh_constants_cuda.h
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-01-21 00:35:49 UTC (rev 19406)
@@ -26,7 +26,7 @@
!=====================================================================
*/
-/* daniel: trivia
+/* trivia
- for most working arrays we use now "realw" instead of "float" type declarations to make it easier to switch
between a real or double precision simulation
@@ -123,8 +123,9 @@
//#define MANUALLY_UNROLLED_LOOPS
// cuda kernel block size for updating displacements/potential (newmark time scheme)
-#define BLOCKSIZE_KERNEL1 128
-#define BLOCKSIZE_KERNEL3 128
+// current hardware: 128 is slightly faster than 256 ( ~ 4%)
+#define BLOCKSIZE_KERNEL1 128
+#define BLOCKSIZE_KERNEL3 128
#define BLOCKSIZE_TRANSFER 256
/* ----------------------------------------------------------------------------------------------- */
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/noise_tomography_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/noise_tomography_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/noise_tomography_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -29,7 +29,11 @@
#include <stdio.h>
#include <cuda.h>
#include <cublas.h>
+
+#ifdef WITH_MPI
#include <mpi.h>
+#endif
+
#include <sys/types.h>
#include <unistd.h>
@@ -56,7 +60,11 @@
TRACE("fortranprint");
int procid;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
+#else
+ procid = 0;
+#endif
printf("%d: sends msg_id %d\n",procid,*id);
}
@@ -67,7 +75,11 @@
TRACE("fortranprintf");
int procid;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
+#else
+ procid = 0;
+#endif
printf("%d: sends val %e\n",procid,*val);
}
@@ -78,7 +90,11 @@
TRACE("fortranprintd");
int procid;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
+#else
+ procid = 0;
+#endif
printf("%d: sends val %e\n",procid,*val);
}
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 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -29,7 +29,10 @@
#include <stdio.h>
#include <cuda.h>
#include <cublas.h>
+
+#ifdef WITH_MPI
#include <mpi.h>
+#endif
#include <sys/time.h>
#include <sys/resource.h>
@@ -67,7 +70,11 @@
void pause_for_debugger(int pause) {
if(pause) {
int myrank;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+#else
+ myrank = 0;
+#endif
printf("I'm rank %d\n",myrank);
int i = 0;
char hostname[256];
@@ -103,7 +110,7 @@
{
printf("\nERROR: %s\n",info);
fflush(stdout);
-#ifdef USE_MPI
+#ifdef WITH_MPI
MPI_Abort(MPI_COMM_WORLD,1);
#endif
//free(info);
@@ -119,7 +126,7 @@
{
printf("\nCUDA error !!!!! <%s> !!!!! \nat CUDA call error code: # %d\n",cudaGetErrorString(err),num);
fflush(stdout);
-#ifdef USE_MPI
+#ifdef WITH_MPI
MPI_Abort(MPI_COMM_WORLD,1);
#endif
exit(EXIT_FAILURE);
@@ -149,9 +156,8 @@
/* ----------------------------------------------------------------------------------------------- */
// Saves GPU memory usage to file
-void output_free_memory(char* info_str) {
- int myrank;
- MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+void output_free_memory(int myrank,char* info_str) {
+
FILE* fp;
char filename[BUFSIZ];
double free_db,used_db,total_db;
@@ -170,21 +176,26 @@
// Fortran-callable version of above method
extern "C"
void FC_FUNC_(output_free_device_memory,
- OUTPUT_FREE_DEVICE_MEMORY)(int* id) {
+ OUTPUT_FREE_DEVICE_MEMORY)(int* myrank) {
TRACE("output_free_device_memory");
char info[6];
- sprintf(info,"f %d:",*id);
- output_free_memory(info);
+ sprintf(info,"f %d:",*myrank);
+ output_free_memory(*myrank,info);
}
/* ----------------------------------------------------------------------------------------------- */
+/*
void show_free_memory(char* info_str) {
// show memory usage of GPU
int myrank;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+#else
+ myrank = 0;
+#endif
double free_db,used_db,total_db;
get_free_memory(&free_db,&used_db,&total_db);
@@ -194,17 +205,18 @@
}
+ extern "C"
+ void FC_FUNC_(show_free_device_memory,
+ SHOW_FREE_DEVICE_MEMORY)() {
+ TRACE("show_free_device_memory");
+
+ show_free_memory("from fortran");
+ }
+ */
+
/* ----------------------------------------------------------------------------------------------- */
-extern "C"
-void FC_FUNC_(show_free_device_memory,
- SHOW_FREE_DEVICE_MEMORY)() {
-TRACE("show_free_device_memory");
- show_free_memory("from fortran");
-}
-
-
extern "C"
void FC_FUNC_(get_free_device_memory,
get_FREE_DEVICE_MEMORY)(realw* free, realw* used, realw* total ) {
@@ -430,7 +442,7 @@
int myrank = *myrank_f;
// cuda initialization (needs -lcuda library)
- // daniel: note, cuInit initializes the driver API.
+ // note: cuInit initializes the driver API.
// it is needed for any following CUDA driver API function call (format cuFUNCTION(..) )
// however, for the CUDA runtime API functions (format cudaFUNCTION(..) )
// the initialization is implicit, thus cuInit() here would not be needed...
@@ -457,7 +469,7 @@
exit_on_error("CUDA Compute capability major number should be at least 1.3\n");
}
- //daniel: from here on we use the runtime API ...
+ // note: from here on we use the runtime API ...
// Gets number of GPU devices
int device_count = 0;
cudaGetDeviceCount(&device_count);
@@ -681,9 +693,6 @@
print_CUDA_error_if_any(cudaMemcpy(mp->d_ispec_is_inner, h_ispec_is_inner,
mp->NSPEC_AB*sizeof(int),cudaMemcpyHostToDevice),1206);
- // daniel: check
- //check_ispec_is(Mesh_pointer,0);
-
// absorbing boundaries
mp->d_num_abs_boundary_faces = *h_num_abs_boundary_faces;
if( *ABSORBING_CONDITIONS && mp->d_num_abs_boundary_faces > 0 ){
@@ -693,10 +702,6 @@
(mp->d_num_abs_boundary_faces)*sizeof(int),
cudaMemcpyHostToDevice),1102);
- // daniel: check
- //check_array_ispec(Mesh_pointer,1);
-
-
print_CUDA_error_if_any(cudaMalloc((void**) &(mp->d_abs_boundary_ijk),
3*NGLL2*(mp->d_num_abs_boundary_faces)*sizeof(int)),1103);
print_CUDA_error_if_any(cudaMemcpy(mp->d_abs_boundary_ijk, h_abs_boundary_ijk,
@@ -1065,9 +1070,6 @@
print_CUDA_error_if_any(cudaMemcpy(mp->d_ispec_is_elastic,ispec_is_elastic,
mp->NSPEC_AB*sizeof(int),cudaMemcpyHostToDevice),4012);
- // daniel: check
- //check_ispec_is(Mesh_pointer_f,1);
-
// phase elements
mp->num_phase_ispec_elastic = *num_phase_ispec_elastic;
print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_phase_ispec_inner_elastic),
@@ -1075,9 +1077,6 @@
print_CUDA_error_if_any(cudaMemcpy(mp->d_phase_ispec_inner_elastic,phase_ispec_inner_elastic,
mp->num_phase_ispec_elastic*2*sizeof(int),cudaMemcpyHostToDevice),4011);
- //daniel: check
- //check_phase_ispec(Mesh_pointer_f,1);
-
// for seismograms
if( mp->nrec_local > 0 ){
print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_station_seismo_field),
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c 2012-01-21 00:35:49 UTC (rev 19406)
@@ -30,14 +30,22 @@
#include <stdlib.h>
#include <math.h>
#include <errno.h>
+
+#ifdef WITH_MPI
#include <mpi.h>
+#endif
+
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
void save_to_max_surface_file_(float* maxval) {
int rank;
char filename[BUFSIZ];
FILE* fp;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+#else
+ rank = 0;
+#endif
sprintf(filename,"maxval_surface_proc_%03d.dat",rank);
fp = fopen(filename,"a+");
fprintf(fp,"%e\n",*maxval);
@@ -80,7 +88,11 @@
int nodes_per_iteration = *nodes_per_iterationf;
char filename[BUFSIZ];
int procid;
+#ifdef WITH_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&procid);
+#else
+ procid = 0;
+#endif
sprintf(filename,"/scratch/eiger/rietmann/SPECFEM3D_AIGLE/in_out_files/DATABASES_MPI/proc%06d_surface_movie",procid);
FILE* fp; int it;
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2012-01-21 00:35:49 UTC (rev 19406)
@@ -1,6 +1,36 @@
-#include "config.h"
+/*
+ !=====================================================================
+ !
+ ! S p e c f e m 3 D V e r s i o n 2 . 0
+ ! ---------------------------------------
+ !
+ ! Main authors: Dimitri Komatitsch and Jeroen Tromp
+ ! Princeton University, USA and University of Pau / CNRS / INRIA
+ ! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+ ! April 2011
+ !
+ ! This program is free software; you can redistribute it and/or modify
+ ! it under the terms of the GNU General Public License as published by
+ ! the Free Software Foundation; either version 2 of the License, or
+ ! (at your option) any later version.
+ !
+ ! This program is distributed in the hope that it will be useful,
+ ! but WITHOUT ANY WARRANTY; without even the implied warranty of
+ ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ ! GNU General Public License for more details.
+ !
+ ! You should have received a copy of the GNU General Public License along
+ ! with this program; if not, write to the Free Software Foundation, Inc.,
+ ! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ !
+ !=====================================================================
+ */
+
#include <stdio.h>
+#include <stdlib.h>
+#include "config.h"
+
typedef float realw;
/* from check_fields_cuda.cu */
@@ -316,7 +346,7 @@
void FC_FUNC_(pause_for_debug,PAUSE_FOR_DEBUG)(){}
void FC_FUNC_(output_free_device_memory,
- OUTPUT_FREE_DEVICE_MEMORY)(int* id){}
+ OUTPUT_FREE_DEVICE_MEMORY)(int* myrank){}
void FC_FUNC_(show_free_device_memory,
SHOW_FREE_DEVICE_MEMORY)(){}
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/write_seismograms_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/write_seismograms_cuda.cu 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/cuda/write_seismograms_cuda.cu 2012-01-21 00:35:49 UTC (rev 19406)
@@ -29,7 +29,7 @@
#include <stdio.h>
#include <cuda.h>
#include <cublas.h>
-#include <mpi.h>
+
#include <sys/types.h>
#include <unistd.h>
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/Makefile.in 2012-01-21 00:35:49 UTC (rev 19406)
@@ -196,10 +196,10 @@
${MPIFCCOMPILE_CHECK} -c -o $O/parallel.o ${SHARED}/parallel.f90
$O/model_external_values.o: ${SHARED}/constants.h model_external_values.f90
- ${MPIFCCOMPILE_CHECK} -c -o $O/model_external_values.o model_external_values.f90
+ ${FCCOMPILE_CHECK} -c -o $O/model_external_values.o model_external_values.f90
$O/model_tomography.o: ${SHARED}/constants.h model_tomography.f90
- ${MPIFCCOMPILE_CHECK} -c -o $O/model_tomography.o model_tomography.f90
+ ${FCCOMPILE_CHECK} -c -o $O/model_tomography.o model_tomography.f90
###
### serial compilation without optimization
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/create_regions_mesh.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -457,7 +457,8 @@
call check_mesh_resolution(myrank,nspec,nglob,ibool,&
xstore_dummy,ystore_dummy,zstore_dummy, &
kappastore,mustore,rho_vp,rho_vs, &
- -1.0d0, model_speed_max,min_resolved_period )
+ -1.0d0, model_speed_max,min_resolved_period, &
+ LOCAL_PATH,SAVE_MESH_FILES )
! saves binary mesh files for attenuation
if( ATTENUATION ) then
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/generate_databases/get_perm_color.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -31,7 +31,7 @@
! start non-blocking MPI calls and overlap them with the calculation of the inner elements
! (which works fine because there are always far more inner elements than outer elements)
-! daniel: modified routines to use element domain flags given in ispec_is_d, thus
+! note: these are modified routines to use element domain flags given in ispec_is_d, thus
! coloring only acoustic or elastic (or..) elements in one run, then repeat run for other domains.
! also, the permutation re-starts at 1 for outer and for inner elements,
! making it usable for the phase_ispec_inner_** arrays for acoustic and elastic elements.
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/check_mesh_resolution.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -27,7 +27,8 @@
subroutine check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
kappastore,mustore,rho_vp,rho_vs, &
- DT, model_speed_max,min_resolved_period )
+ DT, model_speed_max,min_resolved_period, &
+ LOCAL_PATH,SAVE_MESH_FILES )
! check the mesh, stability and resolved period
!
@@ -44,6 +45,9 @@
double precision :: DT
real(kind=CUSTOM_REAL) :: model_speed_max,min_resolved_period
+ character(len=256):: LOCAL_PATH
+ logical :: SAVE_MESH_FILES
+
! local parameters
real(kind=CUSTOM_REAL) :: vpmin,vpmax,vsmin,vsmax,vpmin_glob,vpmax_glob,vsmin_glob,vsmax_glob
real(kind=CUSTOM_REAL) :: distance_min,distance_max,distance_min_glob,distance_max_glob
@@ -57,7 +61,7 @@
integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum
integer :: NGLOB_AB_global_min,NGLOB_AB_global_max,NGLOB_AB_global_sum
integer :: ispec,sizeprocs
-
+
!********************************************************************************
! empirical choice for distorted elements to estimate time step and period resolved:
@@ -73,6 +77,11 @@
logical :: has_vs_zero
real(kind=CUSTOM_REAL),dimension(1) :: tmp_val
+ ! debug: for vtk output
+ real(kind=CUSTOM_REAL),dimension(:),allocatable :: tmp1,tmp2
+ integer:: ier
+ character(len=256) :: filename,prname
+
! initializations
if( DT <= 0.0d0) then
DT_PRESENT = .false.
@@ -99,6 +108,15 @@
has_vs_zero = .false.
+ ! debug: for vtk output
+ if( SAVE_MESH_FILES ) then
+ allocate(tmp1(NSPEC_AB),stat=ier)
+ allocate(tmp2(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp'
+ tmp1(:) = 0.0
+ tmp2(:) = 0.0
+ endif
+
! checks courant number & minimum resolved period for each grid cell
do ispec=1,NSPEC_AB
@@ -133,8 +151,12 @@
if( DT_PRESENT ) then
cmax = max( vpmax,vsmax ) * DT / distance_min
cmax_glob = max(cmax_glob,cmax)
+
+ ! debug: for vtk output
+ if( SAVE_MESH_FILES ) tmp1(ispec) = cmax
endif
-
+
+
! suggested timestep
dt_suggested = COURANT_SUGGESTED * distance_min / max( vpmax,vsmax )
dt_suggested_glob = min( dt_suggested_glob, dt_suggested)
@@ -164,6 +186,9 @@
!pmax = distance_max / min( vpmin,vsmin ) * NELEM_PER_WAVELENGTH
!pmax_glob = max(pmax_glob,pmax)
+ ! debug: for vtk output
+ if( SAVE_MESH_FILES ) tmp2(ispec) = pmax
+
enddo
! determines global min/max values from all cpu partitions
@@ -307,6 +332,24 @@
call bcast_all_cr(tmp_val,1)
min_resolved_period = tmp_val(1)
+ ! debug: for vtk output
+ if( SAVE_MESH_FILES ) then
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ ! courant number
+ if( DT_PRESENT ) then
+ filename = trim(prname)//'res_courant_number'
+ call write_VTK_data_elem_cr(NSPEC_AB,NGLOB_AB, &
+ xstore,ystore,zstore,ibool, &
+ tmp1,filename)
+ endif
+ ! minimum period estimate
+ filename = trim(prname)//'res_minimum_period'
+ call write_VTK_data_elem_cr(NSPEC_AB,NGLOB_AB, &
+ xstore,ystore,zstore,ibool, &
+ tmp2,filename)
+ deallocate(tmp1,tmp2)
+ endif
+
end subroutine check_mesh_resolution
!
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/compute_arrays_source.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/compute_arrays_source.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -321,7 +321,7 @@
!!! ! reads in adjoint source trace
!!! do itime = 1, NSTEP
!!!
-!!! ! things become a bit tricky because of the Newark time scheme at
+!!! ! things become a bit tricky because of the Newmark time scheme at
!!! ! the very beginning of the time loop. however, when we read in the backward/reconstructed
!!! ! wavefields at the end of the first time loop, we can use the adjoint source index from 1 to NSTEP
!!! ! (and then access it in reverse NSTEP-it+1 down to 1, for it=1,..NSTEP; see compute_add_sources*.f90).
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 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/read_parameter_file.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -223,7 +223,8 @@
implicit none
include "constants.h"
- logical GPU_MODE,GRAVITY
+ logical :: GPU_MODE
+ logical :: GRAVITY
! initializes flags
GPU_MODE = .false.
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/serial.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -38,7 +38,16 @@
!
double precision function wtime()
- wtime = 0.d0
+
+ implicit none
+ real :: ct
+
+ ! note: for simplicity, we take cpu_time which returns the elapsed CPU time in seconds
+ ! (instead of wall clock time for parallel MPI function)
+ call cpu_time(ct)
+
+ wtime = ct
+
end function wtime
!
@@ -660,6 +669,39 @@
!----
!
+ subroutine send_dp(sendbuf, sendcount, dest, sendtag)
+
+ implicit none
+
+ integer dest,sendtag
+ integer sendcount
+ double precision,dimension(sendcount):: sendbuf
+
+ stop 'send_dp not implemented for serial code'
+
+ end subroutine send_dp
+
+!
+!----
+!
+
+ subroutine recv_dp(recvbuf, recvcount, dest, recvtag)
+
+ implicit none
+
+ integer dest,recvtag
+ integer recvcount
+ double precision,dimension(recvcount):: recvbuf
+
+ stop 'recv_dp not implemented for serial code'
+
+ end subroutine recv_dp
+
+!
+!----
+!
+
+
subroutine sendv_cr(sendbuf, sendcount, dest, sendtag)
implicit none
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/write_VTK_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/write_VTK_data.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/shared/write_VTK_data.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -478,3 +478,75 @@
end subroutine write_VTK_data_elem_vectors
+!=============================================================
+
+ subroutine write_VTK_data_elem_cr(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ elem_flag,prname_file)
+
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nspec,nglob
+
+ ! global coordinates
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy
+
+ ! element flag array
+ real(kind=CUSTOM_REAL), dimension(nspec) :: elem_flag
+
+ ! file name
+ character(len=256) :: prname_file
+ !character(len=2), optional, intent(in) :: str_id
+
+ ! local parameters
+ integer :: ispec,i
+
+! write source and receiver VTK files for Paraview
+ !debug
+ !write(IMAIN,*) ' vtk file: '
+ !write(IMAIN,*) ' ',prname_file(1:len_trim(prname_file))//'.vtk'
+
+ open(IOVTK,file=prname_file(1:len_trim(prname_file))//'.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 3.1'
+ write(IOVTK,'(a)') 'material model VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+ do i=1,nglob
+ write(IOVTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i)
+ enddo
+ write(IOVTK,*) ""
+
+ ! note: indices for vtk start at 0
+ write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ do ispec=1,nspec
+ write(IOVTK,'(9i12)') 8,ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1,ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1,&
+ ibool(1,1,NGLLZ,ispec)-1,ibool(NGLLX,1,NGLLZ,ispec)-1,ibool(NGLLX,NGLLY,NGLLZ,ispec)-1,ibool(1,NGLLY,NGLLZ,ispec)-1
+ enddo
+ write(IOVTK,*) ""
+
+ ! type: hexahedrons
+ write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec
+ write(IOVTK,*) (12,ispec=1,nspec)
+ write(IOVTK,*) ""
+
+ write(IOVTK,'(a,i12)') "CELL_DATA ",nspec
+ !if( present( str_id ) ) then
+ ! write(IOVTK,'(a)') "SCALARS elem_flag_"//str_id//" integer"
+ !else
+ ! write(IOVTK,'(a)') "SCALARS elem_flag integer"
+ !endif
+ write(IOVTK,'(a)') "SCALARS elem_val float"
+ write(IOVTK,'(a)') "LOOKUP_TABLE default"
+ do ispec = 1,nspec
+ write(IOVTK,*) elem_flag(ispec)
+ enddo
+ write(IOVTK,*) ""
+ close(IOVTK)
+
+
+ end subroutine write_VTK_data_elem_cr
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/Makefile.in 2012-01-21 00:35:49 UTC (rev 19406)
@@ -51,8 +51,8 @@
@COND_CUDA_TRUE at NVCC = nvcc
@COND_CUDA_FALSE at NVCC = g++
- at COND_CUDA_TRUE@NVCC_FLAGS = $(CUDA_INC) $(MPI_INC) -DCUDA -gencode=arch=compute_20,code=sm_20
- at COND_CUDA_FALSE@NVCC_FLAGS = $(MPI_INC)
+ at COND_CUDA_TRUE@NVCC_FLAGS = $(CUDA_INC) $(MPI_INC) $(COND_MPI_CPPFLAGS) -DCUDA -gencode=arch=compute_20,code=sm_20
+ at COND_CUDA_FALSE@NVCC_FLAGS = $(MPI_INC) $(COND_MPI_CPPFLAGS)
FC = @FC@
@@ -252,8 +252,8 @@
xcombine_surf_data: $O/combine_surf_data.shared.o $O/write_c_binary.cc.o $O/param_reader.cc.o
${FCCOMPILE_CHECK} -o ${E}/xcombine_surf_data $O/combine_surf_data.shared.o $O/write_c_binary.cc.o $O/param_reader.cc.o
-xsmooth_vol_data: $O/smooth_vol_data.o $O/read_parameter_file.shared.o $O/read_value_parameters.shared.o $O/get_value_parameters.shared.o $O/param_reader.cc.o $O/gll_library.shared.o $O/exit_mpi.shared.o $O/parallel.o
- ${FCLINK} -o ${E}/xsmooth_vol_data $O/smooth_vol_data.o $O/read_parameter_file.shared.o $O/read_value_parameters.shared.o $O/get_value_parameters.shared.o $O/param_reader.cc.o $O/gll_library.shared.o $O/exit_mpi.shared.o $O/parallel.o $(MPILIBS)
+xsmooth_vol_data: $O/smooth_vol_data.o $O/read_parameter_file.shared.o $O/read_value_parameters.shared.o $O/get_value_parameters.shared.o $O/param_reader.cc.o $O/gll_library.shared.o $O/exit_mpi.shared.o $(COND_MPI_OBJECTS)
+ ${FCLINK} -o ${E}/xsmooth_vol_data $O/smooth_vol_data.o $O/read_parameter_file.shared.o $O/read_value_parameters.shared.o $O/get_value_parameters.shared.o $O/param_reader.cc.o $O/gll_library.shared.o $O/exit_mpi.shared.o $(COND_MPI_OBJECTS) $(MPILIBS)
clean:
@@ -311,6 +311,9 @@
$O/parallel.o: $(SHARED)constants.h ${SHARED}/parallel.f90
${MPIFCCOMPILE_CHECK} -c -o $O/parallel.o ${SHARED}/parallel.f90
+$O/serial.o: $(SHARED)constants.h ${SHARED}/serial.f90
+ ${FCCOMPILE_CHECK} -c -o $O/serial.o ${SHARED}/serial.f90
+
$O/smooth_vol_data.o: $(SHARED)constants.h ${SHARED}/smooth_vol_data.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/smooth_vol_data.o ${SHARED}/smooth_vol_data.f90
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_acoustic.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -235,12 +235,12 @@
!
! backward/reconstructed wavefields:
! time for b_potential( it ) would correspond to (NSTEP - it - 1 )*DT - t0
-! if we read in saved wavefields b_potential() before Newark time scheme
+! if we read in saved wavefields b_potential() before Newmark time scheme
! (see sources for simulation_type 1 and seismograms)
-! since at the beginning of the time loop, the numerical Newark time scheme updates
+! since at the beginning of the time loop, the numerical Newmark time scheme updates
! the wavefields, that is b_potential( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
!
-! b_potential is now read in after Newark time scheme:
+! b_potential is now read in after Newmark time scheme:
! we read the backward/reconstructed wavefield at the end of the first time loop,
! such that b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! assuming that until that end the backward/reconstructed wavefield and adjoint fields
@@ -387,7 +387,7 @@
endif ! nadj_rec_local > 0
endif
-! note: b_potential() is read in after Newark time scheme, thus
+! note: b_potential() is read in after Newmark time scheme, thus
! b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! thus indexing is NSTEP - it , instead of NSTEP - it - 1
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_add_sources_elastic.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -213,12 +213,12 @@
!
! backward/reconstructed wavefields:
! time for b_displ( it ) would correspond to (NSTEP - it - 1 )*DT - t0
-! if we read in saved wavefields b_displ() before Newark time scheme
+! if we read in saved wavefields b_displ() before Newmark time scheme
! (see sources for simulation_type 1 and seismograms)
-! since at the beginning of the time loop, the numerical Newark time scheme updates
+! since at the beginning of the time loop, the numerical Newmark time scheme updates
! the wavefields, that is b_displ( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
!
-! b_displ is now read in after Newark time scheme:
+! b_displ is now read in after Newmark time scheme:
! we read the backward/reconstructed wavefield at the end of the first time loop,
! such that b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! assuming that until that end the backward/reconstructed wavefield and adjoint fields
@@ -370,7 +370,7 @@
endif ! nadj_rec_local
endif !adjoint
-! note: b_displ() is read in after Newark time scheme, thus
+! note: b_displ() is read in after Newmark time scheme, thus
! b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! thus indexing is NSTEP - it , instead of NSTEP - it - 1
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_acoustic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_acoustic_el.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_acoustic_el.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -106,7 +106,7 @@
! continuity of pressure and normal displacement on global point
!
- ! note: newark time scheme together with definition of scalar potential:
+ ! note: Newmark time scheme together with definition of scalar potential:
! pressure = - chi_dot_dot
! requires that this coupling term uses the updated displacement at time step [t+delta_t],
! which is done at the very beginning of the time loop
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_ac.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_coupling_elastic_ac.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -101,7 +101,7 @@
! continuity of displacement and pressure on global point
!
- ! note: newark time scheme together with definition of scalar potential:
+ ! note: Newmark time scheme together with definition of scalar potential:
! pressure = - chi_dot_dot
! requires that this coupling term uses the *UPDATED* pressure (chi_dot_dot), i.e.
! pressure at time step [t + delta_t]
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_acoustic.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -386,7 +386,7 @@
num_PML_ispec,PML_ispec,&
chi1_dot_dot,chi2_t_dot_dot,chi3_dot_dot,chi4_dot_dot)
- ! Newark time scheme corrector terms
+ ! Newmark time scheme corrector terms
call PML_acoustic_time_corrector(NSPEC_AB,ispec_is_acoustic,deltatover2,&
num_PML_ispec,PML_ispec,PML_damping_d,&
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
@@ -395,7 +395,7 @@
! update velocity
-! note: Newark finite-difference time scheme with acoustic domains:
+! note: Newmark finite-difference time scheme with acoustic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_forces_elastic.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -308,7 +308,7 @@
endif
! updates velocities
-! Newark finite-difference time scheme with elastic domains:
+! Newmark finite-difference time scheme with elastic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_acoustic.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -84,7 +84,7 @@
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
! reads in absorbing boundary array when first phase is running
if( phase_is_inner .eqv. .false. ) then
- ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
! uses fortran routine
!read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
!if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_elastic.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/compute_stacey_elastic.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -91,7 +91,7 @@
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
! reads in absorbing boundary array when first phase is running
if( phase_is_inner .eqv. .false. ) then
- ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
! uses fortran routine
!read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
!if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -149,15 +149,18 @@
gammaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for databases'
+
! mesh node locations
allocate(xstore(NGLOB_AB), &
ystore(NGLOB_AB), &
zstore(NGLOB_AB),stat=ier)
- if( ier /= 0 ) stop 'error allocating arrays for mesh nodes'
+ if( ier /= 0 ) stop 'error allocating arrays for mesh nodes'
+
! material properties
allocate(kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for material properties'
+
! material flags
allocate(ispec_is_acoustic(NSPEC_AB), &
ispec_is_elastic(NSPEC_AB), &
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/iterate_time.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -70,7 +70,7 @@
call it_check_stability()
endif
- ! update displacement using Newark time scheme
+ ! update displacement using Newmark time scheme
call it_update_displacement_scheme()
! acoustic solver
@@ -84,7 +84,7 @@
if( POROELASTIC_SIMULATION ) stop 'poroelastic simulation not implemented yet'
! restores last time snapshot saved for backward/reconstruction of wavefields
- ! note: this must be read in after the Newark time scheme
+ ! note: this must be read in after the Newmark time scheme
! note 2: GPU b_field transfers are included
if( SIMULATION_TYPE == 3 .and. it == 1 ) then
call it_read_foward_arrays()
@@ -136,7 +136,7 @@
end subroutine iterate_time
- !=====================================================================
+!=====================================================================
subroutine it_print_elapsed_time()
use specfem_par
@@ -162,6 +162,8 @@
endif
end subroutine it_print_elapsed_time
+!=====================================================================
+
subroutine it_check_stability()
! computes the maximum of the norm of the displacement
@@ -265,7 +267,7 @@
iseconds = int_tCPU - 3600*ihours - 60*iminutes
write(IMAIN,*) 'Elapsed time in seconds = ',tCPU
write(IMAIN,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
- write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
+ write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',sngl(tCPU/dble(it))
if( ELASTIC_SIMULATION ) then
write(IMAIN,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
else
@@ -285,7 +287,7 @@
iseconds_remain = int_t_remain - 3600*ihours_remain - 60*iminutes_remain
write(IMAIN,*) 'Time steps done = ',it,' out of ',NSTEP
write(IMAIN,*) 'Time steps remaining = ',NSTEP - it
- write(IMAIN,*) 'Estimated remaining time in seconds = ',t_remain
+ write(IMAIN,*) 'Estimated remaining time in seconds = ',sngl(t_remain)
write(IMAIN,"(' Estimated remaining time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
ihours_remain,iminutes_remain,iseconds_remain
@@ -295,7 +297,7 @@
ihours_total = int_t_total / 3600
iminutes_total = (int_t_total - 3600*ihours_total) / 60
iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
- write(IMAIN,*) 'Estimated total run time in seconds = ',t_total
+ write(IMAIN,*) 'Estimated total run time in seconds = ',sngl(t_total)
write(IMAIN,"(' Estimated total run time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
ihours_total,iminutes_total,iseconds_total
write(IMAIN,*) 'We have done ',sngl(100.d0*dble(it)/dble(NSTEP)),'% of that'
@@ -342,7 +344,7 @@
subroutine it_update_displacement_scheme()
-! explicit Newark time scheme with acoustic & elastic domains:
+! explicit Newmark time scheme with acoustic & elastic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
@@ -555,7 +557,7 @@
! note backward/reconstructed wavefields:
! storing wavefield displ() at time step it, corresponds to time (it-1)*DT - t0 (see routine write_seismograms_to_file )
! reconstucted wavefield b_displ() at it corresponds to time (NSTEP-it-1)*DT - t0
-! we read in the reconstructed wavefield at the end of the time iteration loop, i.e. after the Newark scheme,
+! we read in the reconstructed wavefield at the end of the time iteration loop, i.e. after the Newmark scheme,
! thus, indexing is NSTEP-it (rather than something like NSTEP-(it-1) )
if (SIMULATION_TYPE == 3 .and. mod(NSTEP-it,NSTEP_Q_SAVE) == 0) then
! reads files content
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 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -102,6 +102,7 @@
! elastic
! number of elastic elements in this partition
nspec_elastic = count(ispec_is_elastic(:))
+
! elastic simulation
call any_all_l( ANY(ispec_is_elastic), ELASTIC_SIMULATION )
if( ELASTIC_SIMULATION ) then
@@ -169,8 +170,9 @@
if( ier /= 0 ) stop 'error allocating array one_minus_sum_beta etc.'
! reads mass matrices
- read(27) rmass
-
+ read(27,iostat=ier) rmass
+ if( ier /= 0 ) stop 'error reading in array rmass'
+
if( OCEANS ) then
! ocean mass matrix
allocate(rmass_ocean_load(NGLOB_AB),stat=ier)
@@ -183,21 +185,32 @@
endif
!pll material parameters for stacey conditions
- read(27) rho_vp
- read(27) rho_vs
-
+ read(27,iostat=ier) rho_vp
+ if( ier /= 0 ) stop 'error reading in array rho_vp'
+ read(27,iostat=ier) rho_vs
+ if( ier /= 0 ) stop 'error reading in array 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
+ rhostore = 0.0_CUSTOM_REAL
where( mustore > TINYVAL )
rhostore = (rho_vs*rho_vs) / mustore
- elsewhere
- rhostore = 0.0_CUSTOM_REAL
endwhere
+
+ ! note: the construct below leads to a segmentation fault (ifort v11.1). not sure why...
+ ! (where statement - standard fortran 95)
+ !where( mustore > TINYVAL )
+ ! rhostore = (rho_vs*rho_vs) / mustore
+ !elsewhere
+ ! rhostore = 0.0_CUSTOM_REAL
+ !endwhere
endif
endif
else
@@ -419,7 +432,8 @@
call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, &
ibool,xstore,ystore,zstore, &
kappastore,mustore,rho_vp,rho_vs, &
- DT,model_speed_max,min_resolved_period )
+ DT,model_speed_max,min_resolved_period, &
+ LOCAL_PATH,SAVE_MESH_FILES )
else if( ACOUSTIC_SIMULATION ) then
allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array rho_vp'
@@ -430,7 +444,8 @@
call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, &
ibool,xstore,ystore,zstore, &
kappastore,mustore,rho_vp,rho_vs, &
- DT,model_speed_max,min_resolved_period )
+ DT,model_speed_max,min_resolved_period, &
+ LOCAL_PATH,SAVE_MESH_FILES )
deallocate(rho_vp,rho_vs)
endif
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/src/specfem3D/specfem3D_par.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -159,8 +159,8 @@
double precision :: DT
logical :: ATTENUATION,USE_OLSEN_ATTENUATION, &
- OCEANS,TOPOGRAPHY,ABSORBING_CONDITIONS,ANISOTROPY, &
- GRAVITY
+ OCEANS,TOPOGRAPHY,ABSORBING_CONDITIONS,ANISOTROPY
+ logical :: GRAVITY
logical :: SAVE_FORWARD,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
Modified: seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/locate_partition.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/locate_partition.f90 2012-01-20 23:41:33 UTC (rev 19405)
+++ seismo/3D/SPECFEM3D/branches/SPECFEM3D_SUNFLOWER/utils/locate_partition.f90 2012-01-21 00:35:49 UTC (rev 19406)
@@ -26,11 +26,15 @@
! utility to locate partition which is closest to given point location
!
-! compile with:
+! compile with one of these (use your default):
! > gfortran -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
+! > ifort -assume byterecl -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
+! > pgf90 -I src/shared -o bin/xlocate_partition utils/locate_partition.f90
!
-! run with:
-! > ../../bin/xlocate_partition -300000.0 11000.0 -3000.0 in_out_files/DATABASES_MPI/
+! specify a target (x,y) may be in UTM, not lon-lat, then run with:
+! > ./bin/xlocate_partition 70000.0 11000.0 -3000.0 ./in_out_files/DATABASES_MPI/
+!
+! this will generate the output file in_out_files/DATABASES_MPI/partition_bounds.dat
program locate_partition
@@ -83,6 +87,9 @@
print *,'----------------------------'
print *
+ ! open a text file to list the maximal bounds of each partition
+ open(11,file=trim(LOCAL_PATH)//'partition_bounds.dat',status='unknown')
+
! loops over slices (process partitions)
total_distance = HUGEVAL
total_partition = -1
@@ -100,20 +107,26 @@
status='old',action='read',form='unformatted',iostat=ios)
if( ios /= 0 ) exit
- read(27) NSPEC_AB
- read(27) NGLOB_AB
+ read(27,iostat=ier) NSPEC_AB
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ read(27,iostat=ier) NGLOB_AB
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
! ibool file
allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibool'
- read(27) ibool
+ read(27,iostat=ier) ibool
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
! global point arrays
allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array xstore etc.'
- read(27) xstore
- read(27) ystore
- read(27) zstore
+ read(27,iostat=ier) xstore
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ read(27,iostat=ier) ystore
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
+ read(27,iostat=ier) zstore
+ if( ier /= 0 ) stop 'please check your compilation, use the same compiler & flags as for SPECFEM3D'
close(27)
print*,'partition: ',iproc
@@ -122,6 +135,8 @@
print*,' min/max z = ',minval(zstore),maxval(zstore)
print*
+ write(11,'(i10,6e18.6)') iproc,minval(xstore),maxval(xstore),minval(ystore),maxval(ystore),minval(zstore),maxval(zstore)
+
! gets distance to target location
call get_closest_point(target_x,target_y,target_z, &
NGLOB_AB,NSPEC_AB,xstore,ystore,zstore,ibool, &
@@ -140,6 +155,8 @@
enddo ! all slices for points
+ close(11)
+
! checks
if (total_partition < 0 ) then
print*,'Error: partition not found among ',iproc,'partitions searched'
More information about the CIG-COMMITS
mailing list