[cig-commits] r20647 - in seismo/3D/SPECFEM3D/trunk/src: cuda decompose_mesh_SCOTCH generate_databases meshfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Aug 29 22:17:39 PDT 2012
Author: danielpeter
Date: 2012-08-29 22:17:38 -0700 (Wed, 29 Aug 2012)
New Revision: 20647
Modified:
seismo/3D/SPECFEM3D/trunk/src/cuda/specfem3D_gpu_cuda_method_stubs.c
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
Log:
updates tomographic model routines and bug fixes decomposition
Modified: seismo/3D/SPECFEM3D/trunk/src/cuda/specfem3D_gpu_cuda_method_stubs.c
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2012-08-30 02:44:18 UTC (rev 20646)
+++ seismo/3D/SPECFEM3D/trunk/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2012-08-30 05:17:38 UTC (rev 20647)
@@ -1,4 +1,4 @@
-/*
+/*
!=====================================================================
!
! S p e c f e m 3 D V e r s i o n 2 . 1
@@ -33,59 +33,59 @@
typedef float realw;
-
+
//
// src/cuda/check_fields_cuda.cu
//
-void FC_FUNC_(pause_for_debug,PAUSE_FOR_DEBUG)() {}
+void FC_FUNC_(pause_for_debug,PAUSE_FOR_DEBUG)() {}
void FC_FUNC_(output_free_device_memory,
- OUTPUT_FREE_DEVICE_MEMORY)(int* myrank) {}
+ OUTPUT_FREE_DEVICE_MEMORY)(int* myrank) {}
void FC_FUNC_(get_free_device_memory,
- get_FREE_DEVICE_MEMORY)(realw* free, realw* used, realw* total ) {}
+ get_FREE_DEVICE_MEMORY)(realw* free, realw* used, realw* total ) {}
void FC_FUNC_(check_max_norm_displ_gpu,
- CHECK_MAX_NORM_DISPL_GPU)(int* size, realw* displ,long* Mesh_pointer_f,int* announceID) {}
+ CHECK_MAX_NORM_DISPL_GPU)(int* size, realw* displ,long* Mesh_pointer_f,int* announceID) {}
void FC_FUNC_(check_max_norm_vector,
- CHECK_MAX_NORM_VECTOR)(int* size, realw* vector1, int* announceID) {}
+ CHECK_MAX_NORM_VECTOR)(int* size, realw* vector1, int* announceID) {}
void FC_FUNC_(check_max_norm_displ,
- CHECK_MAX_NORM_DISPL)(int* size, realw* displ, int* announceID) {}
+ CHECK_MAX_NORM_DISPL)(int* size, realw* displ, int* announceID) {}
void FC_FUNC_(check_max_norm_b_displ_gpu,
- CHECK_MAX_NORM_B_DISPL_GPU)(int* size, realw* b_displ,long* Mesh_pointer_f,int* announceID) {}
+ CHECK_MAX_NORM_B_DISPL_GPU)(int* size, realw* b_displ,long* Mesh_pointer_f,int* announceID) {}
void FC_FUNC_(check_max_norm_b_accel_gpu,
- CHECK_MAX_NORM_B_ACCEL_GPU)(int* size, realw* b_accel,long* Mesh_pointer_f,int* announceID) {}
+ CHECK_MAX_NORM_B_ACCEL_GPU)(int* size, realw* b_accel,long* Mesh_pointer_f,int* announceID) {}
void FC_FUNC_(check_max_norm_b_veloc_gpu,
- CHECK_MAX_NORM_B_VELOC_GPU)(int* size, realw* b_veloc,long* Mesh_pointer_f,int* announceID) {}
+ CHECK_MAX_NORM_B_VELOC_GPU)(int* size, realw* b_veloc,long* Mesh_pointer_f,int* announceID) {}
void FC_FUNC_(check_max_norm_b_displ,
- CHECK_MAX_NORM_B_DISPL)(int* size, realw* b_displ,int* announceID) {}
+ CHECK_MAX_NORM_B_DISPL)(int* size, realw* b_displ,int* announceID) {}
void FC_FUNC_(check_max_norm_b_accel,
- CHECK_MAX_NORM_B_ACCEL)(int* size, realw* b_accel,int* announceID) {}
+ CHECK_MAX_NORM_B_ACCEL)(int* size, realw* b_accel,int* announceID) {}
void FC_FUNC_(check_error_vectors,
- CHECK_ERROR_VECTORS)(int* sizef, realw* vector1,realw* vector2) {}
+ CHECK_ERROR_VECTORS)(int* sizef, realw* vector1,realw* vector2) {}
void FC_FUNC_(get_max_accel,
- GET_MAX_ACCEL)(int* itf,int* sizef,long* Mesh_pointer) {}
+ GET_MAX_ACCEL)(int* itf,int* sizef,long* Mesh_pointer) {}
void FC_FUNC_(get_norm_acoustic_from_device,
GET_NORM_ACOUSTIC_FROM_DEVICE)(realw* norm,
long* Mesh_pointer_f,
- int* SIMULATION_TYPE) {}
+ int* SIMULATION_TYPE) {}
void FC_FUNC_(get_norm_elastic_from_device,
GET_NORM_ELASTIC_FROM_DEVICE)(realw* norm,
long* Mesh_pointer_f,
- int* SIMULATION_TYPE) {}
+ int* SIMULATION_TYPE) {}
//
@@ -98,7 +98,7 @@
int* NSOURCESf,
int* SIMULATION_TYPEf,
double* h_stf_pre_compute,
- int* myrankf) {}
+ int* myrankf) {}
void FC_FUNC_(compute_add_sources_ac_s3_cuda,
COMPUTE_ADD_SOURCES_AC_s3_CUDA)(long* Mesh_pointer_f,
@@ -106,7 +106,7 @@
int* NSOURCESf,
int* SIMULATION_TYPEf,
double* h_stf_pre_compute,
- int* myrankf) {}
+ int* myrankf) {}
void FC_FUNC_(add_sources_ac_sim_2_or_3_cuda,
ADD_SOURCES_AC_SIM_2_OR_3_CUDA)(long* Mesh_pointer,
@@ -120,7 +120,7 @@
int* time_index,
int* h_islice_selected_rec,
int* nadj_rec_local,
- int* NTSTEP_BETWEEN_READ_ADJSRC) {}
+ int* NTSTEP_BETWEEN_READ_ADJSRC) {}
//
@@ -132,21 +132,21 @@
int* phase_is_innerf,
int* NSOURCESf,
double* h_stf_pre_compute,
- int* myrankf) {}
+ int* myrankf) {}
void FC_FUNC_(compute_add_sources_el_s3_cuda,
COMPUTE_ADD_SOURCES_EL_S3_CUDA)(long* Mesh_pointer,
double* h_stf_pre_compute,
int* NSOURCESf,
int* phase_is_inner,
- int* myrank) {}
+ int* myrank) {}
void FC_FUNC_(add_source_master_rec_noise_cu,
ADD_SOURCE_MASTER_REC_NOISE_CU)(long* Mesh_pointer_f,
int* myrank_f,
int* it_f,
int* irec_master_noise_f,
- int* islice_selected_rec) {}
+ int* islice_selected_rec) {}
void FC_FUNC_(add_sources_el_sim_type_2_or_3,
ADD_SOURCES_EL_SIM_TYPE_2_OR_3)(long* Mesh_pointer,
@@ -160,7 +160,7 @@
int* time_index,
int* h_islice_selected_rec,
int* nadj_rec_local,
- int* NTSTEP_BETWEEN_READ_ADJSRC) {}
+ int* NTSTEP_BETWEEN_READ_ADJSRC) {}
//
@@ -171,13 +171,13 @@
COMPUTE_COUPLING_AC_EL_CUDA)(long* Mesh_pointer_f,
int* phase_is_innerf,
int* num_coupling_ac_el_facesf,
- int* SIMULATION_TYPEf) {}
+ int* SIMULATION_TYPEf) {}
void FC_FUNC_(compute_coupling_el_ac_cuda,
COMPUTE_COUPLING_EL_AC_CUDA)(long* Mesh_pointer_f,
int* phase_is_innerf,
int* num_coupling_ac_el_facesf,
- int* SIMULATION_TYPEf) {}
+ int* SIMULATION_TYPEf) {}
//
@@ -194,7 +194,7 @@
int* max_nibool_interfaces_ext_mesh,
int* nibool_interfaces_ext_mesh,
int* ibool_interfaces_ext_mesh,
- int* FORWARD_OR_ADJOINT){}
+ int* FORWARD_OR_ADJOINT){}
void FC_FUNC_(transfer_asmbl_pot_to_device,
TRANSFER_ASMBL_POT_TO_DEVICE)(
@@ -205,31 +205,31 @@
int* max_nibool_interfaces_ext_mesh,
int* nibool_interfaces_ext_mesh,
int* ibool_interfaces_ext_mesh,
- int* FORWARD_OR_ADJOINT) {}
+ int* FORWARD_OR_ADJOINT) {}
void FC_FUNC_(compute_forces_acoustic_cuda,
COMPUTE_FORCES_ACOUSTIC_CUDA)(long* Mesh_pointer_f,
int* iphase,
int* nspec_outer_acoustic,
int* nspec_inner_acoustic,
- int* SIMULATION_TYPE) {}
+ int* SIMULATION_TYPE) {}
void FC_FUNC_(kernel_3_a_acoustic_cuda,KERNEL_3_ACOUSTIC_CUDA)(
long* Mesh_pointer,
int* size_F,
- int* SIMULATION_TYPE) {}
+ int* SIMULATION_TYPE) {}
void FC_FUNC_(kernel_3_b_acoustic_cuda,KERNEL_3_ACOUSTIC_CUDA)(
long* Mesh_pointer,
int* size_F,
realw* deltatover2_F,
int* SIMULATION_TYPE,
- realw* b_deltatover2_F) {}
+ realw* b_deltatover2_F) {}
void FC_FUNC_(acoustic_enforce_free_surf_cuda,
ACOUSTIC_ENFORCE_FREE_SURF_CUDA)(long* Mesh_pointer_f,
int* SIMULATION_TYPE,
- int* ABSORB_FREE_SURFACE) {}
+ int* ABSORB_FREE_SURFACE) {}
//
@@ -243,17 +243,17 @@
int* max_nibool_interfaces_ext_mesh,
int* nibool_interfaces_ext_mesh,
int* ibool_interfaces_ext_mesh,
- int* FORWARD_OR_ADJOINT){}
+ int* FORWARD_OR_ADJOINT){}
void FC_FUNC_(transfer_boundary_from_device_a,
TRANSFER_BOUNDARY_FROM_DEVICE_A)(long* Mesh_pointer,
- int* nspec_outer_elastic) {}
+ int* nspec_outer_elastic) {}
void FC_FUNC_(transfer_boundary_to_device_a,
TRANSFER_BOUNDARY_TO_DEVICE_A)(long* Mesh_pointer,
realw* buffer_recv_vector_ext_mesh,
int* num_interfaces_ext_mesh,
- int* max_nibool_interfaces_ext_mesh) {}
+ int* max_nibool_interfaces_ext_mesh) {}
//void FC_FUNC_(assemble_accel_on_device,
// ASSEMBLE_ACCEL_on_DEVICE)(long* Mesh_pointer, realw* accel,
@@ -262,7 +262,7 @@
// int* max_nibool_interfaces_ext_mesh,
// int* nibool_interfaces_ext_mesh,
// int* ibool_interfaces_ext_mesh,
-// int* FORWARD_OR_ADJOINT) {}
+// int* FORWARD_OR_ADJOINT) {}
void FC_FUNC_(transfer_asmbl_accel_to_device,
TRANSFER_ASMBL_ACCEL_TO_DEVICE)(long* Mesh_pointer, realw* accel,
@@ -271,7 +271,7 @@
int* max_nibool_interfaces_ext_mesh,
int* nibool_interfaces_ext_mesh,
int* ibool_interfaces_ext_mesh,
- int* FORWARD_OR_ADJOINT) {}
+ int* FORWARD_OR_ADJOINT) {}
void FC_FUNC_(compute_forces_elastic_cuda,
COMPUTE_FORCES_ELASTIC_CUDA)(long* Mesh_pointer_f,
@@ -281,12 +281,12 @@
int* SIMULATION_TYPE,
int* COMPUTE_AND_STORE_STRAIN,
int* ATTENUATION,
- int* ANISOTROPY) {}
+ int* ANISOTROPY) {}
void FC_FUNC_(sync_copy_from_device,
SYNC_copy_FROM_DEVICE)(long* Mesh_pointer_f,
int* iphase,
- realw* send_buffer) {}
+ realw* send_buffer) {}
void FC_FUNC_(kernel_3_a_cuda,
KERNEL_3_A_CUDA)(long* Mesh_pointer,
@@ -294,18 +294,18 @@
realw* deltatover2_F,
int* SIMULATION_TYPE_f,
realw* b_deltatover2_F,
- int* OCEANS) {}
+ int* OCEANS) {}
void FC_FUNC_(kernel_3_b_cuda,
KERNEL_3_B_CUDA)(long* Mesh_pointer,
int* size_F,
realw* deltatover2_F,
int* SIMULATION_TYPE_f,
- realw* b_deltatover2_F) {}
+ realw* b_deltatover2_F) {}
void FC_FUNC_(elastic_ocean_load_cuda,
ELASTIC_OCEAN_LOAD_CUDA)(long* Mesh_pointer_f,
- int* SIMULATION_TYPE) {}
+ int* SIMULATION_TYPE) {}
//
@@ -314,23 +314,23 @@
void FC_FUNC_(compute_kernels_elastic_cuda,
COMPUTE_KERNELS_ELASTIC_CUDA)(long* Mesh_pointer,
- realw* deltat_f) {}
+ realw* deltat_f) {}
void FC_FUNC_(compute_kernels_strgth_noise_cu,
COMPUTE_KERNELS_STRGTH_NOISE_CU)(long* Mesh_pointer,
realw* h_noise_surface_movie,
- realw* deltat) {}
+ realw* deltat) {}
void FC_FUNC_(compute_kernels_acoustic_cuda,
COMPUTE_KERNELS_ACOUSTIC_CUDA)(
long* Mesh_pointer,
- realw* deltat_f) {}
+ realw* deltat_f) {}
void FC_FUNC_(compute_kernels_hess_cuda,
COMPUTE_KERNELS_HESS_CUDA)(long* Mesh_pointer,
realw* deltat_f,
int* ELASTIC_SIMULATION,
- int* ACOUSTIC_SIMULATION) {}
+ int* ACOUSTIC_SIMULATION) {}
//
@@ -343,7 +343,7 @@
int* phase_is_innerf,
int* SIMULATION_TYPEf,
int* SAVE_FORWARDf,
- realw* h_b_absorb_potential) {}
+ realw* h_b_absorb_potential) {}
//
@@ -355,7 +355,7 @@
int* phase_is_innerf,
int* SIMULATION_TYPEf,
int* SAVE_FORWARDf,
- realw* h_b_absorb_field) {}
+ realw* h_b_absorb_field) {}
//
@@ -363,10 +363,10 @@
//
void FC_FUNC_(initialize_cuda_device,
- INITIALIZE_CUDA_DEVICE)(int* myrank_f,int* ncuda_devices) {
+ INITIALIZE_CUDA_DEVICE)(int* myrank_f,int* ncuda_devices) {
fprintf(stderr,"ERROR: GPU_MODE enabled without GPU/CUDA Support. To enable GPU support, reconfigure with --with-cuda flag.\n");
exit(1);
-}
+}
//
@@ -382,7 +382,7 @@
int* SIMULATION_TYPE,
realw* b_deltat_F,
realw* b_deltatsqover2_F,
- realw* b_deltatover2_F) {}
+ realw* b_deltatover2_F) {}
void FC_FUNC_(it_update_displacement_ac_cuda,
it_update_displacement_ac_cuda)(long* Mesh_pointer_f,
@@ -393,31 +393,31 @@
int* SIMULATION_TYPE,
realw* b_deltat_F,
realw* b_deltatsqover2_F,
- realw* b_deltatover2_F) {}
+ realw* b_deltatover2_F) {}
//
// src/cuda/noise_tomography_cuda.cu
//
-void FC_FUNC_(fortranflush,FORTRANFLUSH)(int* rank){}
+void FC_FUNC_(fortranflush,FORTRANFLUSH)(int* rank){}
-void FC_FUNC_(fortranprint,FORTRANPRINT)(int* id) {}
+void FC_FUNC_(fortranprint,FORTRANPRINT)(int* id) {}
-void FC_FUNC_(fortranprintf,FORTRANPRINTF)(realw* val) {}
+void FC_FUNC_(fortranprintf,FORTRANPRINTF)(realw* val) {}
-void FC_FUNC_(fortranprintd,FORTRANPRINTD)(double* val) {}
+void FC_FUNC_(fortranprintd,FORTRANPRINTD)(double* val) {}
-void FC_FUNC_(make_displ_rand,MAKE_DISPL_RAND)(long* Mesh_pointer_f,realw* h_displ) {}
+void FC_FUNC_(make_displ_rand,MAKE_DISPL_RAND)(long* Mesh_pointer_f,realw* h_displ) {}
void FC_FUNC_(transfer_surface_to_host,
TRANSFER_SURFACE_TO_HOST)(long* Mesh_pointer_f,
- realw* h_noise_surface_movie) {}
+ realw* h_noise_surface_movie) {}
void FC_FUNC_(noise_read_add_surface_movie_cu,
NOISE_READ_ADD_SURFACE_MOVIE_CU)(long* Mesh_pointer_f,
realw* h_noise_surface_movie,
- int* NOISE_TOMOGRAPHYf) {}
+ int* NOISE_TOMOGRAPHYf) {}
//
@@ -461,7 +461,7 @@
int* request_send_vector_ext_mesh,
int* request_recv_vector_ext_mesh,
realw* buffer_recv_vector_ext_mesh
- ) {}
+ ) {}
void FC_FUNC_(prepare_fields_acoustic_device,
PREPARE_FIELDS_ACOUSTIC_DEVICE)(long* Mesh_pointer_f,
@@ -486,12 +486,12 @@
realw* coupling_ac_el_jacobian2Dw,
int* num_colors_outer_acoustic,
int* num_colors_inner_acoustic,
- int* num_elem_colors_acoustic) {}
+ int* num_elem_colors_acoustic) {}
void FC_FUNC_(prepare_fields_acoustic_adj_dev,
PREPARE_FIELDS_ACOUSTIC_ADJ_DEV)(long* Mesh_pointer_f,
int* SIMULATION_TYPE,
- int* APPROXIMATE_HESS_KL) {}
+ int* APPROXIMATE_HESS_KL) {}
void FC_FUNC_(prepare_fields_elastic_device,
PREPARE_FIELDS_ELASTIC_DEVICE)(long* Mesh_pointer_f,
@@ -546,7 +546,7 @@
realw *c46store,
realw *c55store,
realw *c56store,
- realw *c66store){}
+ realw *c66store){}
void FC_FUNC_(prepare_fields_elastic_adj_dev,
PREPARE_FIELDS_ELASTIC_ADJ_DEV)(long* Mesh_pointer_f,
@@ -561,7 +561,7 @@
int* R_size,
realw* b_R_xx,realw* b_R_yy,realw* b_R_xy,realw* b_R_xz,realw* b_R_yz,
realw* b_alphaval,realw* b_betaval,realw* b_gammaval,
- int* APPROXIMATE_HESS_KL){}
+ int* APPROXIMATE_HESS_KL){}
void FC_FUNC_(prepare_sim2_or_3_const_device,
PREPARE_SIM2_OR_3_CONST_DEVICE)(
@@ -570,7 +570,7 @@
int* islice_selected_rec_size,
int* nadj_rec_local,
int* nrec,
- int* myrank) {}
+ int* myrank) {}
void FC_FUNC_(prepare_fields_noise_device,
PREPARE_FIELDS_NOISE_DEVICE)(long* Mesh_pointer_f,
@@ -586,7 +586,7 @@
realw* normal_y_noise,
realw* normal_z_noise,
realw* mask_noise,
- realw* free_surface_jacobian2Dw) {}
+ realw* free_surface_jacobian2Dw) {}
void FC_FUNC_(prepare_fields_gravity_device,
PREPARE_FIELDS_gravity_DEVICE)(long* Mesh_pointer_f,
@@ -595,10 +595,10 @@
realw* minus_g,
realw* h_wgll_cube,
int* ACOUSTIC_SIMULATION,
- realw* rhostore) {}
+ realw* rhostore) {}
void FC_FUNC_(prepare_seismogram_fields,
- PREPARE_SEISMOGRAM_FIELDS)(long* Mesh_pointer,int* nrec_local, double* nu, double* hxir, double* hetar, double* hgammar) {}
+ PREPARE_SEISMOGRAM_FIELDS)(long* Mesh_pointer,int* nrec_local, double* nu, double* hxir, double* hetar, double* hgammar) {}
void FC_FUNC_(prepare_cleanup_device,
PREPARE_CLEANUP_DEVICE)(long* Mesh_pointer_f,
@@ -612,7 +612,7 @@
int* ATTENUATION,
int* ANISOTROPY,
int* OCEANS,
- int* APPROXIMATE_HESS_KL) {}
+ int* APPROXIMATE_HESS_KL) {}
//
@@ -620,41 +620,41 @@
//
void FC_FUNC_(transfer_fields_el_to_device,
- TRANSFER_FIELDS_EL_TO_DEVICE)(int* size, realw* displ, realw* veloc, realw* accel,long* Mesh_pointer_f) {}
+ TRANSFER_FIELDS_EL_TO_DEVICE)(int* size, realw* displ, realw* veloc, realw* accel,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_fields_el_from_device,
- TRANSFER_FIELDS_EL_FROM_DEVICE)(int* size, realw* displ, realw* veloc, realw* accel,long* Mesh_pointer_f) {}
+ TRANSFER_FIELDS_EL_FROM_DEVICE)(int* size, realw* displ, realw* veloc, realw* accel,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_fields_to_device,
TRANSFER_B_FIELDS_TO_DEVICE)(int* size, realw* b_displ, realw* b_veloc, realw* b_accel,
- long* Mesh_pointer_f) {}
+ long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_fields_from_device,
- TRANSFER_B_FIELDS_FROM_DEVICE)(int* size, realw* b_displ, realw* b_veloc, realw* b_accel,long* Mesh_pointer_f) {}
+ TRANSFER_B_FIELDS_FROM_DEVICE)(int* size, realw* b_displ, realw* b_veloc, realw* b_accel,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_accel_to_device,
- TRNASFER_ACCEL_TO_DEVICE)(int* size, realw* accel,long* Mesh_pointer_f) {}
+ TRNASFER_ACCEL_TO_DEVICE)(int* size, realw* accel,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_accel_from_device,
- TRANSFER_ACCEL_FROM_DEVICE)(int* size, realw* accel,long* Mesh_pointer_f) {}
+ TRANSFER_ACCEL_FROM_DEVICE)(int* size, realw* accel,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_accel_from_device,
- TRNASFER_B_ACCEL_FROM_DEVICE)(int* size, realw* b_accel,long* Mesh_pointer_f) {}
+ TRNASFER_B_ACCEL_FROM_DEVICE)(int* size, realw* b_accel,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_sigma_from_device,
- TRANSFER_SIGMA_FROM_DEVICE)(int* size, realw* sigma_kl,long* Mesh_pointer_f) {}
+ TRANSFER_SIGMA_FROM_DEVICE)(int* size, realw* sigma_kl,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_displ_from_device,
- TRANSFER_B_DISPL_FROM_DEVICE)(int* size, realw* displ,long* Mesh_pointer_f) {}
+ TRANSFER_B_DISPL_FROM_DEVICE)(int* size, realw* displ,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_displ_from_device,
- TRANSFER_DISPL_FROM_DEVICE)(int* size, realw* displ,long* Mesh_pointer_f) {}
+ TRANSFER_DISPL_FROM_DEVICE)(int* size, realw* displ,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_compute_kernel_answers_from_device,
TRANSFER_COMPUTE_KERNEL_ANSWERS_FROM_DEVICE)(long* Mesh_pointer,
realw* rho_kl,int* size_rho,
realw* mu_kl, int* size_mu,
- realw* kappa_kl, int* size_kappa) {}
+ realw* kappa_kl, int* size_kappa) {}
void FC_FUNC_(transfer_compute_kernel_fields_from_device,
TRANSFER_COMPUTE_KERNEL_FIELDS_FROM_DEVICE)(long* Mesh_pointer,
@@ -677,7 +677,7 @@
realw* kappa_kl, int* size_kappa,
realw* epsilon_trace_over_3,
realw* b_epsilon_trace_over_3,
- int* size_epsilon_trace_over_3) {}
+ int* size_epsilon_trace_over_3) {}
void FC_FUNC_(transfer_b_fields_att_to_device,
TRANSFER_B_FIELDS_ATT_TO_DEVICE)(long* Mesh_pointer,
@@ -688,7 +688,7 @@
realw* b_epsilondev_xy,
realw* b_epsilondev_xz,
realw* b_epsilondev_yz,
- int* size_epsilondev) {}
+ int* size_epsilondev) {}
void FC_FUNC_(transfer_fields_att_from_device,
TRANSFER_FIELDS_ATT_FROM_DEVICE)(long* Mesh_pointer,
@@ -699,19 +699,19 @@
realw* epsilondev_xy,
realw* epsilondev_xz,
realw* epsilondev_yz,
- int* size_epsilondev) {}
+ int* size_epsilondev) {}
void FC_FUNC_(transfer_kernels_el_to_host,
TRANSFER_KERNELS_EL_TO_HOST)(long* Mesh_pointer,
realw* h_rho_kl,
realw* h_mu_kl,
realw* h_kappa_kl,
- int* NSPEC_AB) {}
+ int* NSPEC_AB) {}
void FC_FUNC_(transfer_kernels_noise_to_host,
TRANSFER_KERNELS_NOISE_TO_HOST)(long* Mesh_pointer,
realw* h_Sigma_kl,
- int* NSPEC_AB) {}
+ int* NSPEC_AB) {}
void FC_FUNC_(transfer_fields_ac_to_device,
TRANSFER_FIELDS_AC_TO_DEVICE)(
@@ -719,7 +719,7 @@
realw* potential_acoustic,
realw* potential_dot_acoustic,
realw* potential_dot_dot_acoustic,
- long* Mesh_pointer_f) {}
+ long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_fields_ac_to_device,
TRANSFER_B_FIELDS_AC_TO_DEVICE)(
@@ -727,14 +727,14 @@
realw* b_potential_acoustic,
realw* b_potential_dot_acoustic,
realw* b_potential_dot_dot_acoustic,
- long* Mesh_pointer_f) {}
+ long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_fields_ac_from_device,
TRANSFER_FIELDS_AC_FROM_DEVICE)(int* size,
realw* potential_acoustic,
realw* potential_dot_acoustic,
realw* potential_dot_dot_acoustic,
- long* Mesh_pointer_f) {}
+ long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_fields_ac_from_device,
TRANSFER_B_FIELDS_AC_FROM_DEVICE)(
@@ -742,29 +742,29 @@
realw* b_potential_acoustic,
realw* b_potential_dot_acoustic,
realw* b_potential_dot_dot_acoustic,
- long* Mesh_pointer_f) {}
+ long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_dot_dot_from_device,
- TRNASFER_DOT_DOT_FROM_DEVICE)(int* size, realw* potential_dot_dot_acoustic,long* Mesh_pointer_f) {}
+ TRNASFER_DOT_DOT_FROM_DEVICE)(int* size, realw* potential_dot_dot_acoustic,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_b_dot_dot_from_device,
- TRNASFER_B_DOT_DOT_FROM_DEVICE)(int* size, realw* b_potential_dot_dot_acoustic,long* Mesh_pointer_f) {}
+ TRNASFER_B_DOT_DOT_FROM_DEVICE)(int* size, realw* b_potential_dot_dot_acoustic,long* Mesh_pointer_f) {}
void FC_FUNC_(transfer_kernels_ac_to_host,
TRANSFER_KERNELS_AC_TO_HOST)(long* Mesh_pointer,
realw* h_rho_ac_kl,
realw* h_kappa_ac_kl,
- int* NSPEC_AB) {}
+ int* NSPEC_AB) {}
void FC_FUNC_(transfer_kernels_hess_el_tohost,
TRANSFER_KERNELS_HESS_EL_TOHOST)(long* Mesh_pointer,
realw* h_hess_kl,
- int* NSPEC_AB) {}
+ int* NSPEC_AB) {}
void FC_FUNC_(transfer_kernels_hess_ac_tohost,
TRANSFER_KERNELS_HESS_AC_TOHOST)(long* Mesh_pointer,
realw* h_hess_ac_kl,
- int* NSPEC_AB) {}
+ int* NSPEC_AB) {}
//
@@ -778,14 +778,14 @@
realw* seismograms_d,
realw* seismograms_v,
realw* seismograms_a,
- int* it) {}
+ int* it) {}
void FC_FUNC_(transfer_station_el_from_device,
TRANSFER_STATION_EL_FROM_DEVICE)(realw* displ,realw* veloc,realw* accel,
realw* b_displ, realw* b_veloc, realw* b_accel,
long* Mesh_pointer_f,int* number_receiver_global,
int* ispec_selected_rec,int* ispec_selected_source,
- int* ibool,int* SIMULATION_TYPEf) {}
+ int* ibool,int* SIMULATION_TYPEf) {}
void FC_FUNC_(transfer_station_ac_from_device,
TRANSFER_STATION_AC_FROM_DEVICE)(
@@ -800,5 +800,5 @@
int* ispec_selected_rec,
int* ispec_selected_source,
int* ibool,
- int* SIMULATION_TYPEf) {}
+ int* SIMULATION_TYPEf) {}
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-08-30 02:44:18 UTC (rev 20646)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-08-30 05:17:38 UTC (rev 20647)
@@ -679,6 +679,8 @@
! gets materials id associations
allocate(num_material(1:nspec),stat=ier)
if( ier /= 0 ) stop 'error allocating array num_material'
+ ! note: num_material can be negative for tomographic material elements
+ ! (which are counted then as elastic elements)
num_material(:) = mat(1,:)
! in case of acoustic/elastic/poro simulation, weights elements accordingly
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-08-30 02:44:18 UTC (rev 20646)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-08-30 05:17:38 UTC (rev 20647)
@@ -1274,18 +1274,24 @@
! sets weights for elements
do el = 0, nspec-1
+ ! note: num_material index can be negative for tomographic material definitions
! acoustic element (cheap)
- if ( is_acoustic(num_material(el+1)) ) then
- elmnts_load(el+1) = elmnts_load(el+1)*ACOUSTIC_LOAD
- endif
- ! elastic element (expensive)
- if ( is_elastic(num_material(el+1)) ) then
+ if( num_material(el+1) > 0 ) then
+ if ( is_acoustic(num_material(el+1)) ) then
+ elmnts_load(el+1) = elmnts_load(el+1)*ACOUSTIC_LOAD
+ endif
+ ! elastic element (expensive)
+ if ( is_elastic(num_material(el+1)) ) then
+ elmnts_load(el+1) = elmnts_load(el+1)*ELASTIC_LOAD
+ endif
+ ! poroelastic element (very expensive)
+ if ( is_poroelastic(num_material(el+1)) ) then
+ elmnts_load(el+1) = elmnts_load(el+1)*POROELASTIC_LOAD
+ endif
+ else
+ ! tomographic materials count as elastic
elmnts_load(el+1) = elmnts_load(el+1)*ELASTIC_LOAD
endif
- ! poroelastic element (very expensive)
- if ( is_poroelastic(num_material(el+1)) ) then
- elmnts_load(el+1) = elmnts_load(el+1)*POROELASTIC_LOAD
- endif
enddo
end subroutine acoustic_elastic_poro_load
@@ -1363,13 +1369,17 @@
! counts coupled elements
nfaces_coupled = 0
do el = 0, nspec-1
- if ( is_poroelastic(num_material(el+1)) ) then
+ if( num_material(el+1) > 0 ) then
+ if ( is_poroelastic(num_material(el+1)) ) then
do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ if(num_material(adjncy(el_adj)+1) > 0 ) then
+ if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
nfaces_coupled = nfaces_coupled + 1
- endif
+ endif
+ endif
enddo
- endif
+ endif
+ endif
enddo
! coupled elements
@@ -1380,15 +1390,19 @@
! stores elements indices
nfaces_coupled = 0
do el = 0, nspec-1
- if ( is_poroelastic(num_material(el+1)) ) then
+ if( num_material(el+1) > 0 ) then
+ if ( is_poroelastic(num_material(el+1)) ) then
do el_adj = xadj(el), xadj(el+1) - 1
- if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ if( num_material(adjncy(el_adj)+1) > 0 ) then
+ if ( is_elastic(abs(num_material(adjncy(el_adj)+1))) ) then
nfaces_coupled = nfaces_coupled + 1
faces_coupled(1,nfaces_coupled) = el
faces_coupled(2,nfaces_coupled) = adjncy(el_adj)
- endif
+ endif
+ endif
enddo
- endif
+ endif
+ endif
enddo
! puts coupled elements into same partition
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90 2012-08-30 02:44:18 UTC (rev 20646)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90 2012-08-30 05:17:38 UTC (rev 20647)
@@ -85,12 +85,10 @@
!debug
!print*,"nundefMat_ext_mesh:",nundefMat_ext_mesh
-! prepares tomography model if needed for elements with undefined material definitions
- ! TODO: Max -- somehow this code is breaking when I try to run
- ! Piero's PREM
- ! if( nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO ) then
- ! call model_tomography_broadcast(myrank)
- ! endif
+ ! prepares tomography model if needed for elements with undefined material definitions
+ if( nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO ) then
+ call model_tomography_broadcast(myrank)
+ endif
! prepares external model values if needed
select case( IMODEL )
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90 2012-08-30 02:44:18 UTC (rev 20646)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90 2012-08-30 05:17:38 UTC (rev 20647)
@@ -34,12 +34,12 @@
!
!--------------------------------------------------------------------------------------------------
- module tomography
+ module tomography_par
include "constants.h"
! for external tomography:
- ! file must be in ../in_data/files/ directory
+ ! file must be in ../in_data_files/ directory
! (regular spaced, xyz-block file in ascii)
character (len=80) :: TOMO_FILENAME = 'tomography_model.xyz'
@@ -59,14 +59,19 @@
! min/max statistics
double precision :: VP_MIN,VS_MIN,RHO_MIN,VP_MAX,VS_MAX,RHO_MAX
- end module tomography
+ ! process rank
+ integer :: myrank_tomo
+ end module tomography_par
+
!
!-------------------------------------------------------------------------------------------------
!
subroutine model_tomography_broadcast(myrank)
+ use tomography_par, only: myrank_tomo
+
implicit none
! include "constants.h"
@@ -89,6 +94,9 @@
!endif
!call bcast_all_cr(vp_tomography,size(vp_tomography))
+ ! stores rank
+ myrank_tomo = myrank
+
end subroutine model_tomography_broadcast
!
@@ -105,7 +113,7 @@
! this could be problematic for example if the tomographic regions have different resolution
! leading to a waste of memory and cpu time in the partitioning process
- use tomography
+ use tomography_par
implicit none
@@ -115,6 +123,7 @@
real(kind=CUSTOM_REAL) :: x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo
integer :: irecord,ier
character(len=256):: filename
+ character(len=256):: string_read
!TOMO_FILENAME='DATA/veryfast_tomography_abruzzo_complete.xyz'
! probably the simple position for the filename is the constat.h
@@ -124,15 +133,34 @@
! it is a possible solution )
! magnoni 1/12/09
filename = IN_DATA_FILES_PATH(1:len_trim(IN_DATA_FILES_PATH))//trim(TOMO_FILENAME)
+
+ ! user output
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' tomography model : ',trim(filename)
+ endif
+
+ ! opens file for reading
open(unit=27,file=trim(filename),status='old',action='read',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error reading tomography file')
! reads in model dimensions
- read(27,*) ORIG_X, ORIG_Y, ORIG_Z, END_X, END_Y, END_Z
- read(27,*) SPACING_X, SPACING_Y, SPACING_Z
- read(27,*) NX, NY, NZ
- read(27,*) VP_MIN, VP_MAX, VS_MIN, VS_MAX, RHO_MIN, RHO_MAX
+ !read(27,*) ORIG_X, ORIG_Y, ORIG_Z, END_X, END_Y, END_Z
+ call tomo_read_next_line(27,string_read)
+ read(string_read,*) ORIG_X, ORIG_Y, ORIG_Z, END_X, END_Y, END_Z
+ !read(27,*) SPACING_X, SPACING_Y, SPACING_Z
+ call tomo_read_next_line(27,string_read)
+ read(string_read,*) SPACING_X, SPACING_Y, SPACING_Z
+
+ !read(27,*) NX, NY, NZ
+ call tomo_read_next_line(27,string_read)
+ read(string_read,*) NX, NY, NZ
+
+ !read(27,*) VP_MIN, VP_MAX, VS_MIN, VS_MAX, RHO_MIN, RHO_MAX
+ call tomo_read_next_line(27,string_read)
+ read(string_read,*) VP_MIN, VP_MAX, VS_MIN, VS_MAX, RHO_MIN, RHO_MAX
+
+ ! total number of element records
nrecord = NX*NY*NZ
! allocates model records
@@ -142,8 +170,19 @@
z_tomography(1:nrecord),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+ ! first record
+ !read(27,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo
+ call tomo_read_next_line(27,string_read)
+ read(string_read,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo
+
+ ! stores record values
+ vp_tomography(1) = vp_tomo
+ vs_tomography(1) = vs_tomo
+ rho_tomography(1) = rho_tomo
+ z_tomography(1) = z_tomo
+
! reads in record sections
- do irecord = 1,nrecord
+ do irecord = 2,nrecord
read(27,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo
! stores record values
@@ -156,11 +195,49 @@
! user output
if( myrank == 0 ) then
- write(IMAIN,*) ' tomography model: ',trim(TOMO_FILENAME)
+ write(IMAIN,*) ' number of records: ',nrecord
endif
end subroutine read_model_tomography
+!--------------------
+
+ subroutine tomo_read_next_line(unit_in,string_read)
+
+ implicit none
+
+ integer :: unit_in
+ character(len=256) :: string_read
+
+ integer :: ios
+
+ do
+ read(unit=unit_in,fmt="(a256)",iostat=ios) string_read
+ if(ios /= 0) stop 'error while reading tomography file'
+
+ ! suppress leading white spaces, if any
+ string_read = adjustl(string_read)
+
+ ! suppress trailing carriage return (ASCII code 13) if any (e.g. if input text file coming from Windows/DOS)
+ if(index(string_read,achar(13)) > 0) string_read = string_read(1:index(string_read,achar(13))-1)
+
+ ! exit loop when we find the first line that is not a comment or a white line
+ if(len_trim(string_read) == 0) cycle
+ if(string_read(1:1) /= '#') exit
+ enddo
+
+ ! suppress trailing white spaces, if any
+ string_read = string_read(1:len_trim(string_read))
+
+ ! suppress trailing comments, if any
+ if(index(string_read,'#') > 0) string_read = string_read(1:index(string_read,'#')-1)
+
+ ! suppress leading and trailing white spaces again, if any, after having suppressed the leading junk
+ string_read = adjustl(string_read)
+ string_read = string_read(1:len_trim(string_read))
+
+ end subroutine tomo_read_next_line
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -168,15 +245,10 @@
subroutine model_tomography(xmesh,ymesh,zmesh,rho_final,vp_final,vs_final,qmu_atten)
- use tomography
+ use tomography_par
implicit none
- !integer, intent(in) :: NX,NY,NZ
- !real(kind=CUSTOM_REAL), dimension(1:NX*NY*NZ), intent(in) :: vp_tomography,vs_tomography,rho_tomography,z_tomography
- !double precision, intent(in) :: ORIG_X,ORIG_Y,ORIG_Z,SPACING_X,SPACING_Y,SPACING_Z
- !double precision, intent(in) :: VP_MIN,VS_MIN,RHO_MIN,VP_MAX,VS_MAX,RHO_MAX
-
double precision, intent(in) :: xmesh,ymesh,zmesh
real(kind=CUSTOM_REAL), intent(out) :: vp_final,vs_final,rho_final,qmu_atten
@@ -243,6 +315,14 @@
p6 = (ix+1)+(iy+1)*NX+(iz+1)*(NX*NY)
p7 = ix+(iy+1)*NX+(iz+1)*(NX*NY)
+ if(p0 < 0 .or. p1 < 0 .or. p2 < 0 .or. p3 < 0 .or. p4 < 0 .or. p5 < 0 .or. p6 < 0 .or. p7 < 0) then
+ print*,'error rank ',myrank_tomo
+ print*,'corner index:',p0,p1,p2,p3,p4,p5,p6,p7
+ print*,'location:',sngl(xmesh),sngl(ymesh),sngl(zmesh)
+ print*,'origin :',sngl(ORIG_X),sngl(ORIG_Y),sngl(ORIG_Z)
+ call exit_MPI(myrank_tomo,'error corner index in tomography routine')
+ endif
+
if(z_tomography(p4+1) == z_tomography(p0+1)) then
gamma_interp_z1 = 1.d0
else
@@ -294,10 +374,10 @@
gamma_interp_z4 = 0.d0
endif
- gamma_interp_z5 = 1. - gamma_interp_z1
- gamma_interp_z6 = 1. - gamma_interp_z2
- gamma_interp_z7 = 1. - gamma_interp_z3
- gamma_interp_z8 = 1. - gamma_interp_z4
+ gamma_interp_z5 = 1.d0 - gamma_interp_z1
+ gamma_interp_z6 = 1.d0 - gamma_interp_z2
+ gamma_interp_z7 = 1.d0 - gamma_interp_z3
+ gamma_interp_z8 = 1.d0 - gamma_interp_z4
vp1 = vp_tomography(p0+1)
vp2 = vp_tomography(p1+1)
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2012-08-30 02:44:18 UTC (rev 20646)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2012-08-30 05:17:38 UTC (rev 20647)
@@ -73,9 +73,9 @@
double precision, dimension(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) :: xstore,ystore,zstore
- integer ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec)
+ integer,dimension(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) :: ibool
- character(len=256) LOCAL_PATH
+ character(len=256) :: LOCAL_PATH
! auxiliary variables to generate the mesh
! real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: nodes_coords
@@ -85,7 +85,7 @@
integer iax,iay,iar
integer isubregion,nsubregions,doubling_index,nmeshregions
integer imaterial_number
- integer true_material_num(nspec)
+ integer, dimension(nspec) :: true_material_num
integer, dimension(:,:,:), allocatable :: material_num
!integer material_num(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
@@ -157,6 +157,8 @@
integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick
double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick
+ logical, parameter :: DEBUG = .true.
+
! **************
! create the name for the database of the current slide and region
@@ -380,8 +382,61 @@
! note: if mesh squeezes elements such that we can't distinguish two close-by mesh points anymore
! the total number of mesh points might have changed
if(nglob /= NGLOB_AB) then
- print*,'error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB
- call exit_MPI(myrank,'incorrect global number, please check mesh input parameters')
+ ! user output
+ print*,'error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB
+ if(myrank == 0 ) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB
+ write(IMAIN,*)
+ write(IMAIN,*) 'writing out problematic mesh: mesh_debug.vtk'
+ write(IMAIN,*) 'please check your mesh setup...'
+ endif
+ ! debug file output
+ ! vtk file format output
+ open(66,file=prname(1:len_trim(prname))//'mesh_debug.vtk',status='unknown')
+ write(66,'(a)') '# vtk DataFile Version 3.1'
+ write(66,'(a)') 'material model VTK file'
+ write(66,'(a)') 'ASCII'
+ write(66,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(66, '(a,i12,a)') 'POINTS ', nspec*NGLLX_M*NGLLY_M*NGLLZ_M, ' float'
+ ilocnum = 0
+ ibool(:,:,:,:) = 0
+ do ispec=1,nspec
+ do k=1,NGLLZ_M
+ do j=1,NGLLY_M
+ do i=1,NGLLX_M
+ ilocnum = ilocnum + 1
+ ibool(i,j,k,ispec) = ilocnum
+ write(66,*) xstore(i,j,k,ispec),ystore(i,j,k,ispec),zstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+ write(66,*)
+ ! note: indices for vtk start at 0
+ write(66,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+ do ispec=1,nspec
+ write(66,'(9i12)') 8, &
+ ibool(1,1,1,ispec)-1,ibool(2,1,1,ispec)-1,ibool(2,2,1,ispec)-1,ibool(1,2,1,ispec)-1,&
+ ibool(1,1,2,ispec)-1,ibool(2,1,2,ispec)-1,ibool(2,2,2,ispec)-1,ibool(1,2,2,ispec)-1
+ enddo
+ ibool(:,:,:,:) = 0
+ write(66,*)
+ ! type: hexahedrons
+ write(66,'(a,i12)') "CELL_TYPES ",nspec
+ write(66,*) (12,ispec=1,nspec)
+ write(66,*)
+ write(66,'(a,i12)') "CELL_DATA ",nspec
+ write(66,'(a)') "SCALARS elem_val float"
+ write(66,'(a)') "LOOKUP_TABLE default"
+ do ispec = 1,nspec
+ write(66,*) true_material_num(ispec)
+ enddo
+ write(66,*) ""
+ close(66)
+
+ ! stop mesher
+ call exit_MPI(myrank,'incorrect global number, please check mesh input parameters')
end if
! put in classical format
More information about the CIG-COMMITS
mailing list