[cig-commits] r12545 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: . OUTPUT_FILES setup src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Aug 6 07:41:43 PDT 2008
Author: dkomati1
Date: 2008-08-06 07:41:41 -0700 (Wed, 06 Aug 2008)
New Revision: 12545
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_perm_cuthill_mckee.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.bgl
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.rang
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declar.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_1D_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_xi.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90
Log:
completely removed Cuthill-McKee sorting, whose implementation was unsafe
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-08-06 14:41:41 UTC (rev 12545)
@@ -114,7 +114,6 @@
$O/get_jacobian_boundaries.o \
$O/get_jacobian_discontinuities.o \
$O/get_model.o \
- $O/get_perm_cuthill_mckee.o \
$O/get_shape2D.o \
$O/get_shape3D.o \
$O/get_value_parameters.o \
@@ -178,7 +177,7 @@
DEFAULT = \
xcreate_header_file \
$(OUTPUT_FILES_INC)/values_from_mesher.h \
- xmeshfem3D \
+ xspecfem3D \
$(EMPTY_MACRO)
default: $(DEFAULT)
@@ -197,9 +196,9 @@
# rules for the main programs
XMESHFEM_OBJECTS = $O/meshfem3D.o $O/exit_mpi.o $(SOLVER_ARRAY_OBJECTS) $(LIBSPECFEM)
-xmeshfem3D: $(XMESHFEM_OBJECTS)
+xspecfem3D: $(XMESHFEM_OBJECTS)
## use MPI here
- ${MPIFCCOMPILE_CHECK} -o $(BIN)/xmeshfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
+ ${MPIFCCOMPILE_CHECK} -o $(BIN)/xspecfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
# solver also depends on values from mesher
XSPECFEM_OBJECTS = $(SOLVER_ARRAY_OBJECTS) $O/exit_mpi.o $(LIBSPECFEM)
@@ -211,7 +210,7 @@
${MPIFCCOMPILE_CHECK} -o $(BIN)/xcreate_header_file $O/create_header_file.o $O/exit_mpi.o $O/get_value_parameters.o $O/read_compute_parameters.o $O/memory_eval.o $O/save_header_file.o $O/count_number_of_sources.o $O/read_value_parameters.o $O/euler_angles.o $O/reduce.o $O/rthetaphi_xyz.o $O/auto_ner.o
clean:
- rm -f $O/* *.o work.pc* *.mod $(BIN)/xmeshfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file PI*
+ rm -f $O/* *.o work.pc* *.mod $(BIN)/xspecfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file PI*
###
@@ -452,9 +451,6 @@
$O/create_name_database.o: $(SPECINC)/constants.h $S/create_name_database.f90
${FCCOMPILE_CHECK} -c -o $O/create_name_database.o ${FCFLAGS_f90} $S/create_name_database.f90
-$O/get_perm_cuthill_mckee.o: $(SPECINC)/constants.h $S/get_perm_cuthill_mckee.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_perm_cuthill_mckee.o ${FCFLAGS_f90} $S/get_perm_cuthill_mckee.f90
-
$O/define_derivation_matrices.o: $(SPECINC)/constants.h $S/define_derivation_matrices.f90
${FCCOMPILE_CHECK} -c -o $O/define_derivation_matrices.o ${FCFLAGS_f90} $S/define_derivation_matrices.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.bgl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.bgl 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.bgl 2008-08-06 14:41:41 UTC (rev 12545)
@@ -157,7 +157,7 @@
# default targets
DEFAULT = \
OUTPUT_FILES/values_from_mesher.h \
- xmeshfem3D \
+ xspecfem3D \
$(EMPTY_MACRO)
default: $(DEFAULT)
@@ -176,9 +176,9 @@
# rules for the main programs
XMESHFEM_OBJECTS = $O/meshfem3D.o $O/exit_mpi.o $(SOLVER_ARRAY_OBJECTS) $(LIBSPECFEM)
-xmeshfem3D: $(XMESHFEM_OBJECTS)
+xspecfem3D: $(XMESHFEM_OBJECTS)
## use MPI here
- ${MPIFCCOMPILE_CHECK} -o $(BIN)/xmeshfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
+ ${MPIFCCOMPILE_CHECK} -o $(BIN)/xspecfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
# solver also depends on values from mesher
XSPECFEM_OBJECTS = $(SOLVER_ARRAY_OBJECTS) $O/exit_mpi.o $(LIBSPECFEM)
@@ -190,7 +190,7 @@
${MPIFCCOMPILE_CHECK} -o $(BIN)/xcreate_header_file $O/create_header_file.o $O/exit_mpi.o $O/get_value_parameters.o $O/read_compute_parameters.o $O/memory_eval.o $O/save_header_file.o $O/count_number_of_sources.o $O/read_value_parameters.o $O/euler_angles.o $O/reduce.o $O/rthetaphi_xyz.o $O/auto_ner.o
clean:
- rm -f $O/* *.o work.pc* *.mod $(BIN)/xmeshfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file PI*
+ rm -f $O/* *.o work.pc* *.mod $(BIN)/xspecfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file PI*
###
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.rang
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.rang 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile.rang 2008-08-06 14:41:41 UTC (rev 12545)
@@ -158,7 +158,7 @@
DEFAULT = \
xcreate_header_file \
OUTPUT_FILES/values_from_mesher.h \
- xmeshfem3D \
+ xspecfem3D \
$(EMPTY_MACRO)
default: $(DEFAULT)
@@ -177,9 +177,9 @@
# rules for the main programs
XMESHFEM_OBJECTS = $O/meshfem3D.o $O/exit_mpi.o $(SOLVER_ARRAY_OBJECTS) $(LIBSPECFEM)
-xmeshfem3D: $(XMESHFEM_OBJECTS)
+xspecfem3D: $(XMESHFEM_OBJECTS)
## use MPI here
- ${MPIFCCOMPILE_CHECK} -o $(BIN)/xmeshfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
+ ${MPIFCCOMPILE_CHECK} -o $(BIN)/xspecfem3D $(XMESHFEM_OBJECTS) $(MPILIBS)
# solver also depends on values from mesher
XSPECFEM_OBJECTS = $(SOLVER_ARRAY_OBJECTS) $O/exit_mpi.o $(LIBSPECFEM)
@@ -191,7 +191,7 @@
${MPIFCCOMPILE_CHECK} -o $(BIN)/xcreate_header_file $O/create_header_file.o $O/exit_mpi.o $O/get_value_parameters.o $O/read_compute_parameters.o $O/memory_eval.o $O/save_header_file.o $O/count_number_of_sources.o $O/read_value_parameters.o $O/euler_angles.o $O/reduce.o $O/rthetaphi_xyz.o $O/auto_ner.o
clean:
- rm -f $O/* *.o work.pc* *.mod $(BIN)/xmeshfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file PI*
+ rm -f $O/* *.o work.pc* *.mod $(BIN)/xspecfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file PI*
###
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h 2008-08-06 14:41:41 UTC (rev 12545)
@@ -1,4 +1,4 @@
-
+
!
! this is the parameter file for static compilation of the solver
!
@@ -24,57 +24,57 @@
! ---------------------------
!
! exact total number of spectral elements in entire mesh =
- ! 224640.000000000
+ ! 224640.00000000000
! approximate total number of points in entire mesh =
- ! 15038376.0000000
+ ! 15038376.000000000
! approximate total number of degrees of freedom in entire mesh =
- ! 42811848.0000000
+ ! 42811848.000000000
!
! position of the mesh chunk at the surface:
! -----------------------------------------
!
- ! angular size in first direction in degrees = 90.00000
- ! angular size in second direction in degrees = 90.00000
+ ! angular size in first direction in degrees = 90.000000
+ ! angular size in second direction in degrees = 90.000000
!
- ! longitude of center in degrees = 0.0000000E+00
- ! latitude of center in degrees = 90.00000
+ ! longitude of center in degrees = 0.0000000
+ ! latitude of center in degrees = 90.000000
!
- ! angle of rotation of the first chunk = 0.0000000E+00
+ ! angle of rotation of the first chunk = 0.0000000
!
! corner 1
- ! longitude in degrees = -45.0000000000000
- ! latitude in degrees = 35.4465752495873
+ ! longitude in degrees = -45.000000000000000
+ ! latitude in degrees = 35.446575249587262
!
! corner 2
- ! longitude in degrees = 45.0000000000000
- ! latitude in degrees = 35.4465752495873
+ ! longitude in degrees = 45.000000000000000
+ ! latitude in degrees = 35.446575249587262
!
! corner 3
- ! longitude in degrees = -135.000000000000
- ! latitude in degrees = 35.4465752495873
+ ! longitude in degrees = -135.00000000000000
+ ! latitude in degrees = 35.446575249587262
!
! corner 4
- ! longitude in degrees = 135.000000000000
- ! latitude in degrees = 35.4465752495873
+ ! longitude in degrees = 135.00000000000000
+ ! latitude in degrees = 35.446575249587262
!
! resolution of the mesh at the surface:
! -------------------------------------
!
! spectral elements along a great circle = 256
! GLL points along a great circle = 1024
- ! average distance between points in degrees = 0.3515625
- ! average distance between points in km = 39.09196
- ! average size of a spectral element in km = 156.3679
+ ! average distance between points in degrees = 0.35156250
+ ! average distance between points in km = 39.091965
+ ! average size of a spectral element in km = 156.36786
!
! number of time steps = 7900
!
! number of seismic sources = 1
!
-
+
! approximate static memory needed by the solver:
! ----------------------------------------------
!
- ! size of static arrays per slice = 8.134651184082031E-002 GB
+ ! size of static arrays per slice = 8.13465118408203125E-002 GB
!
! (should be below and typically equal to 80% of 1.5 GB = 1.2 GB on pangu
! at Caltech, and below and typically equal to 85% of 2 GB = 1.7 GB
@@ -82,39 +82,39 @@
! (if significantly more, the job will not run by lack of memory)
! (if significantly less, you waste a significant amount of memory)
!
- ! size of static arrays for all slices = 0.325386047363281 GB
- ! = 3.177598118782043E-004 TB
+ ! size of static arrays for all slices = 0.32538604736328125 GB
+ ! = 3.17759811878204346E-004 TB
!
-
+
integer, parameter :: NEX_XI_VAL = 64
integer, parameter :: NEX_ETA_VAL = 64
-
+
integer, parameter :: NSPEC_CRUST_MANTLE = 8640
integer, parameter :: NSPEC_OUTER_CORE = 688
integer, parameter :: NSPEC_INNER_CORE = 32
-
+
integer, parameter :: NGLOB_CRUST_MANTLE = 576013
integer, parameter :: NGLOB_OUTER_CORE = 47985
integer, parameter :: NGLOB_INNER_CORE = 2601
-
+
integer, parameter :: NSPECMAX_ANISO_IC = 1
-
+
integer, parameter :: NSPECMAX_ISO_MANTLE = 8640
integer, parameter :: NSPECMAX_TISO_MANTLE = 1
integer, parameter :: NSPECMAX_ANISO_MANTLE = 1
-
+
integer, parameter :: NSPEC_CRUST_MANTLE_ATTENUAT = 1
integer, parameter :: NSPEC_INNER_CORE_ATTENUATION = 1
-
+
integer, parameter :: NSPEC_CRUST_MANTLE_STR_OR_ATT = 1
integer, parameter :: NSPEC_INNER_CORE_STR_OR_ATT = 1
-
+
integer, parameter :: NSPEC_CRUST_MANTLE_STR_AND_ATT = 1
integer, parameter :: NSPEC_INNER_CORE_STR_AND_ATT = 1
-
+
integer, parameter :: NSPEC_CRUST_MANTLE_STRAIN_ONLY = 1
integer, parameter :: NSPEC_INNER_CORE_STRAIN_ONLY = 1
-
+
integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT = 1
integer, parameter :: NSPEC_OUTER_CORE_ADJOINT = 1
integer, parameter :: NSPEC_INNER_CORE_ADJOINT = 1
@@ -122,29 +122,29 @@
integer, parameter :: NGLOB_OUTER_CORE_ADJOINT = 1
integer, parameter :: NGLOB_INNER_CORE_ADJOINT = 1
integer, parameter :: NSPEC_OUTER_CORE_ROT_ADJOINT = 1
-
+
integer, parameter :: NSPEC_CRUST_MANTLE_STACEY = 1
integer, parameter :: NSPEC_OUTER_CORE_STACEY = 1
-
+
integer, parameter :: NGLOB_CRUST_MANTLE_OCEANS = 1
-
+
logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .false.
-
+
logical, parameter :: ANISOTROPIC_3D_MANTLE_VAL = .false.
-
+
logical, parameter :: ANISOTROPIC_INNER_CORE_VAL = .false.
-
+
logical, parameter :: ATTENUATION_VAL = .false.
-
+
logical, parameter :: ATTENUATION_3D_VAL = .false.
-
+
logical, parameter :: ELLIPTICITY_VAL = .false.
-
+
logical, parameter :: GRAVITY_VAL = .false.
-
+
logical, parameter :: ROTATION_VAL = .false.
integer, parameter :: NSPEC_OUTER_CORE_ROTATION = 1
-
+
integer, parameter :: NGLOB1D_RADIAL_CM = 109
integer, parameter :: NGLOB1D_RADIAL_OC = 65
integer, parameter :: NGLOB1D_RADIAL_IC = 9
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h 2008-08-06 14:41:41 UTC (rev 12545)
@@ -479,11 +479,3 @@
integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
-! for Cuthill-McKee (1969) permutation
- logical, parameter :: PERFORM_CUTHILL_MCKEE = .false.
- integer, parameter :: NGNOD_HEXAHEDRA = 8
-! perform classical or multi-level Cuthill-McKee ordering
- logical, parameter :: CMcK_MULTI = .false.
-! maximum size if multi-level Cuthill-McKee ordering
- integer, parameter :: LIMIT_MULTI_CUTHILL = 50
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -18,13 +18,7 @@
weight = wxgll(i)*wygll(j)*wzgll(k)
-!! DK DK changed this for merged version
-! if(PERFORM_CUTHILL_MCKEE) then
-! iglobnum = ibool(i,j,k,invperm(ispec))
-! iglobnum = ibool(i,j,k,perm(ispec))
-! else
- iglobnum = ibool(i,j,k,ispec)
-! endif
+ iglobnum = ibool(i,j,k,ispec)
! compute the jacobian
xixl = xixstore(i,j,k)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -34,7 +34,6 @@
ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE,NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
- NGLOB2DMAX_XY, &
myrank,LOCAL_PATH,OCEANS,ibathy_topo,rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD,&
ATTENUATION,ATTENUATION_3D, &
NCHUNKS,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
@@ -60,7 +59,7 @@
normal_xmin,normal_xmax,normal_ymin, &
normal_ymax,normal_bottom,normal_top, &
kappavstore,kappahstore,muvstore,muhstore,eta_anisostore,rmass,xelm_store,yelm_store,zelm_store, &
- npoin2D_xi,npoin2D_eta,perm,invperm)
+ npoin2D_xi,npoin2D_eta)
! create the different regions of the mesh
@@ -526,8 +525,6 @@
integer NUMBER_OF_MESH_LAYERS,layer_shift,cpt,first_layer_aniso,last_layer_aniso,FIRST_ELT_NON_ANISO
double precision, dimension(:,:), allocatable :: stretch_tab
- integer :: NGLOB2DMAX_XY
-
integer :: nb_layer_above_aniso,FIRST_ELT_ABOVE_ANISO
integer, parameter :: maxker=200
@@ -577,20 +574,6 @@
double precision r_moho,r_400,r_670
logical :: is_superbrick
-! added for Cuthill McKee permutation
-! integer, dimension(:), allocatable :: perm,perm_tmp,temp_array_1D_int
-!! DK DK added this for merged version, as a quick patch
- integer, dimension(nspec) :: perm,invperm
- integer, dimension(:), allocatable :: perm_tmp,temp_array_1D_int
- logical, dimension(:,:), allocatable :: temp_array_2D_log
- integer, dimension(:,:,:,:), allocatable :: temp_array_int
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: temp_array_real
- double precision, dimension(:,:,:,:), allocatable :: temp_array_dble
- double precision, dimension(:,:,:,:,:), allocatable :: temp_array_dble_5dim
-!! DK DK added this for merged version
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: temp_array_xelm_yelm_zelm
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: temp_rmass
-
! the height at which the central cube is cut
integer :: nz_inf_limit
@@ -1851,11 +1834,11 @@
! create MPI buffers
! arrays locval(npointot) and ifseg(npointot) used to save memory
call get_MPI_cutplanes_xi(myrank,nspec,iMPIcut_xi,ibool, &
- xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_ETA_FACE,iregion_code,NGLOB2DMAX_XY,nglob,iboolleft_xi,iboolright_xi,NGLOB2DMAX_XMIN_XMAX,npoin2D_xi)
+ ifseg,npointot, &
+ NSPEC2D_ETA_FACE,iregion_code,nglob,iboolleft_xi,iboolright_xi,NGLOB2DMAX_XMIN_XMAX,npoin2D_xi)
call get_MPI_cutplanes_eta(myrank,nspec,iMPIcut_eta,ibool, &
- xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_XI_FACE,iregion_code,NGLOB2DMAX_XY,nglob,iboolleft_eta,iboolright_eta,NGLOB2DMAX_YMIN_YMAX,npoin2D_eta)
+ ifseg,npointot, &
+ NSPEC2D_XI_FACE,iregion_code,nglob,iboolleft_eta,iboolright_eta,NGLOB2DMAX_YMIN_YMAX,npoin2D_eta)
call get_MPI_1D_buffers(myrank,nspec,iMPIcut_xi,iMPIcut_eta,ibool,idoubling, &
xstore,ystore,zstore,ifseg,npointot, &
NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,iregion_code,nglob, &
@@ -1876,63 +1859,6 @@
print *,"ERROR can not deallocate ifseg in create_regions_mesh ier=",ier
endif
-!! DK DK for merged code, copied here to be able to create mass matrix in right order
- if (PERFORM_CUTHILL_MCKEE) then
-!!!!!!!!!!!!!!!!!!!! allocate(perm(nspec))
- if(iregion_code == IREGION_CRUST_MANTLE) then
- ! do not permute anisotropic elements
- perm(1:FIRST_ELT_NON_ANISO-1) = (/ (i,i=1,FIRST_ELT_NON_ANISO-1) /)
-
- ! no more connectivity between layers below and above the anisotropic layers => 2 permutations
- ! permute the bottom of the region : below the aniso layers
- allocate(perm_tmp(FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate perm_tmp in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
- call get_perm(ibool(:,:,:,FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1),perm_tmp,LIMIT_MULTI_CUTHILL,&
- (FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO),nglob,.true.,.false.)
- perm(FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1) = perm_tmp(:)+(FIRST_ELT_NON_ANISO-1)
- deallocate(perm_tmp,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate perm_tmp in create_regions_mesh ier=",ier
- endif
-
- ! permute the top of the region : above the aniso layers
- allocate(perm_tmp(nspec-FIRST_ELT_ABOVE_ANISO+1),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate perm_tmp in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- call get_perm(ibool(:,:,:,FIRST_ELT_ABOVE_ANISO:nspec),perm_tmp,LIMIT_MULTI_CUTHILL,&
- (nspec-FIRST_ELT_ABOVE_ANISO+1),nglob,.true.,.false.)
- perm(FIRST_ELT_ABOVE_ANISO:nspec) = perm_tmp(:)+(FIRST_ELT_ABOVE_ANISO-1)
- deallocate(perm_tmp,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate perm_tmp in create_regions_mesh ier=",ier
- endif
- else
- ! the 3 last parameters are : PERFORM_CUTHILL_MCKEE,INVERSE,FACE
- call get_perm(ibool,perm,LIMIT_MULTI_CUTHILL,nspec,nglob,.true.,.false.)
- endif
-!!!!!!!!!!!!!!!!!!!! deallocate(perm)
-!! DK DK create inverse permutation
- if(minval(perm) /= 1) stop 'minval of perm must be 1'
- if(maxval(perm) /= nspec) stop 'maxval of perm must be nspec'
- invperm(:) = -1
- do ispec=1,nspec
- if(invperm(perm(ispec)) == -1) then
- invperm(perm(ispec)) = ispec
- else
- stop 'value already found, permutation is not bijective'
- endif
- enddo
- if(minval(invperm) /= 1) stop 'minval of invperm must be 1'
- if(maxval(invperm) /= nspec) stop 'maxval of invperm must be nspec'
- endif
-!! DK DK for merged code, copied here to be able to create mass matrix in right order
-
! only create mass matrix and save all the final arrays in the second pass
else if(ipass == 2) then
@@ -1947,266 +1873,6 @@
nspec_tiso = 0
endif
-! ***************************************************
-! Cuthill McKee permutation
-! ***************************************************
- if (PERFORM_CUTHILL_MCKEE) then
-!!!!!!!!!!!!!!!!!!! allocate(perm(nspec))
-! if(iregion_code == IREGION_CRUST_MANTLE) then
-! ! do not permute anisotropic elements
-! perm(1:FIRST_ELT_NON_ANISO-1) = (/ (i,i=1,FIRST_ELT_NON_ANISO-1) /)
-
-! ! no more connectivity between layers below and above the anisotropic layers => 2 permutations
-! ! permute the bottom of the region : below the aniso layers
-! allocate(perm_tmp(FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO))
-! call get_perm(ibool(:,:,:,FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1),perm_tmp,LIMIT_MULTI_CUTHILL,&
-! (FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO),nglob,.true.,.false.)
-! perm(FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1) = perm_tmp(:)+(FIRST_ELT_NON_ANISO-1)
-! deallocate(perm_tmp)
-
-! ! permute the top of the region : above the aniso layers
-! allocate(perm_tmp(nspec-FIRST_ELT_ABOVE_ANISO+1))
-! call get_perm(ibool(:,:,:,FIRST_ELT_ABOVE_ANISO:nspec),perm_tmp,LIMIT_MULTI_CUTHILL,&
-! (nspec-FIRST_ELT_ABOVE_ANISO+1),nglob,.true.,.false.)
-! perm(FIRST_ELT_ABOVE_ANISO:nspec) = perm_tmp(:)+(FIRST_ELT_ABOVE_ANISO-1)
-! deallocate(perm_tmp)
-! else
-! ! the 3 last parameters are : PERFORM_CUTHILL_MCKEE,INVERSE,FACE
-! call get_perm(ibool,perm,LIMIT_MULTI_CUTHILL,nspec,nglob,.true.,.false.)
-! endif
-
- ! permutation of xstore, ystore, zstore, rhostore, kappavstore, kappahstore,
- ! muvstore, muhstore, eta_anisostore
-
- allocate(temp_array_dble(NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_dble in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- if(ATTENUATION .and. ATTENUATION_3D) then
- call permute_elements_dble(Qmu_store,temp_array_dble,perm,nspec)
- allocate(temp_array_dble_5dim(N_SLS,NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate permute_elements_dble in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- temp_array_dble_5dim(:,:,:,:,:) = tau_e_store(:,:,:,:,:)
- do i = 1,nspec
- tau_e_store(:,:,:,:,perm(i)) = temp_array_dble_5dim(:,:,:,:,i)
- enddo
- deallocate(temp_array_dble_5dim,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_dble_5dim in create_regions_mesh ier=",ier
- endif
-
- endif
- call permute_elements_dble(xstore,temp_array_dble,perm,nspec)
- call permute_elements_dble(ystore,temp_array_dble,perm,nspec)
- call permute_elements_dble(zstore,temp_array_dble,perm,nspec)
- deallocate(temp_array_dble,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_dble in create_regions_mesh ier=",ier
- endif
-
-!! DK DK added this for merged code
- allocate(temp_array_xelm_yelm_zelm(NGNOD,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_xelm_yelm_zelm in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- call permute_elements_xelm_yelm_zelm(xelm_store,temp_array_xelm_yelm_zelm,perm,nspec)
- call permute_elements_xelm_yelm_zelm(yelm_store,temp_array_xelm_yelm_zelm,perm,nspec)
- call permute_elements_xelm_yelm_zelm(zelm_store,temp_array_xelm_yelm_zelm,perm,nspec)
-
- deallocate(temp_array_xelm_yelm_zelm,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_xelm_yelm_zelm in create_regions_mesh ier=",ier
- endif
-
- allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_real in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- if(ABSORBING_CONDITIONS .and. NCHUNKS /= 6) then
- call permute_elements_real(rho_vp,temp_array_real,perm,nspec)
- call permute_elements_real(rho_vs,temp_array_real,perm,nspec)
- endif
- if((ANISOTROPIC_INNER_CORE .and. iregion_code == IREGION_INNER_CORE) .or. &
- (ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE)) then
- call permute_elements_real(c11store,temp_array_real,perm,nspec)
- call permute_elements_real(c12store,temp_array_real,perm,nspec)
- call permute_elements_real(c13store,temp_array_real,perm,nspec)
- call permute_elements_real(c14store,temp_array_real,perm,nspec)
- call permute_elements_real(c15store,temp_array_real,perm,nspec)
- call permute_elements_real(c16store,temp_array_real,perm,nspec)
- call permute_elements_real(c22store,temp_array_real,perm,nspec)
- call permute_elements_real(c23store,temp_array_real,perm,nspec)
- call permute_elements_real(c24store,temp_array_real,perm,nspec)
- call permute_elements_real(c25store,temp_array_real,perm,nspec)
- call permute_elements_real(c26store,temp_array_real,perm,nspec)
- call permute_elements_real(c33store,temp_array_real,perm,nspec)
- call permute_elements_real(c34store,temp_array_real,perm,nspec)
- call permute_elements_real(c35store,temp_array_real,perm,nspec)
- call permute_elements_real(c36store,temp_array_real,perm,nspec)
- call permute_elements_real(c44store,temp_array_real,perm,nspec)
- call permute_elements_real(c45store,temp_array_real,perm,nspec)
- call permute_elements_real(c46store,temp_array_real,perm,nspec)
- call permute_elements_real(c55store,temp_array_real,perm,nspec)
- call permute_elements_real(c56store,temp_array_real,perm,nspec)
- call permute_elements_real(c66store,temp_array_real,perm,nspec)
- endif
-!! DK DK added this for merged version
- if(iregion_code /= IREGION_OUTER_CORE) then
- call permute_elements_real(kappavstore,temp_array_real,perm,nspec)
- call permute_elements_real(muvstore,temp_array_real,perm,nspec)
- if(iregion_code == IREGION_CRUST_MANTLE .and. NSPECMAX_TISO_MANTLE > 1) then
- if(minval(perm(1:NSPECMAX_TISO_MANTLE)) /= 1) stop 'minval perm for aniso should be 1'
- if(maxval(perm(1:NSPECMAX_TISO_MANTLE)) /= NSPECMAX_TISO_MANTLE) &
- stop 'maxval perm for aniso should be NSPECMAX_TISO_MANTLE'
- call permute_elements_real(kappahstore,temp_array_real,perm,NSPECMAX_TISO_MANTLE)
- call permute_elements_real(muhstore,temp_array_real,perm,NSPECMAX_TISO_MANTLE)
- call permute_elements_real(eta_anisostore,temp_array_real,perm,NSPECMAX_TISO_MANTLE)
- endif
- endif
-
-!! DK DK added this for merged version: attempt to permute mass matrix
-! 333333333333333
- allocate(temp_rmass(NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_rmass in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-! allocate(temp_rmass(nglob))
- do ispec = 1,nspec
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
-! temp_rmass(ibool(i,j,k,perm(ispec))) = rmass(ibool(i,j,k,ispec))
- temp_rmass(i,j,k,ispec) = rmass(ibool(i,j,k,ispec))
- enddo
- enddo
- enddo
- enddo
-! rmass(:) = temp_rmass(:)
- call permute_elements_real(temp_rmass,temp_array_real,perm,nspec)
-! do ispec = 1,nspec
-! do k = 1,NGLLZ
-! do j = 1,NGLLY
-! do i = 1,NGLLX
-! rmass(ibool(i,j,k,ispec)) = temp_rmass(i,j,k,ispec)
-! enddo
-! enddo
-! enddo
-! enddo
-! deallocate(temp_rmass)
-
- deallocate(temp_array_real,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_real in create_regions_mesh ier=",ier
- endif
-
- ! permutation of ibool
- allocate(temp_array_int(NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_int in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- call permute_elements_integer(ibool,temp_array_int,perm,nspec)
- deallocate(temp_array_int,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_int in create_regions_mesh ier=",ier
- endif
-
-!! DK DK added this for merged version: attempt to permute mass matrix
-! 333333333333333
- do ispec = 1,nspec
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- rmass(ibool(i,j,k,ispec)) = temp_rmass(i,j,k,ispec)
- enddo
- enddo
- enddo
- enddo
- deallocate(temp_rmass,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_rmass in create_regions_mesh ier=",ier
- endif
-
-! deallocate(temp_array_real)
-
- ! permutation of iMPIcut_*
- allocate(temp_array_2D_log(2,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_2D_log in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- temp_array_2D_log(:,:) = iMPIcut_xi(:,:)
- do i = 1,nspec
- iMPIcut_xi(:,perm(i)) = temp_array_2D_log(:,i)
- enddo
- temp_array_2D_log(:,:) = iMPIcut_eta(:,:)
- do i = 1,nspec
- iMPIcut_eta(:,perm(i)) = temp_array_2D_log(:,i)
- enddo
-
- deallocate(temp_array_2D_log,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_2D_log in create_regions_mesh ier=",ier
- endif
-
- ! permutation of iboun
- allocate(temp_array_2D_log(6,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_2D_log in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- temp_array_2D_log(:,:) = iboun(:,:)
- do i = 1,nspec
- iboun(:,perm(i)) = temp_array_2D_log(:,i)
- enddo
-
- deallocate(temp_array_2D_log,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_2D_log in create_regions_mesh ier=",ier
- endif
-
-
- ! permutation of idoubling
-!! DK DK added this for merged version because array not allocated in the outer core
- if(iregion_code /= IREGION_OUTER_CORE) then
- allocate(temp_array_1D_int(nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate temp_array_1D_int in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- temp_array_1D_int(:) = idoubling(:)
- do i = 1,nspec
- idoubling(perm(i)) = temp_array_1D_int(i)
- enddo
- deallocate(temp_array_1D_int,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate temp_array_1D_int in create_regions_mesh ier=",ier
- endif
- endif
-
-!!!!!!!!!!!!!!!!!!! deallocate(perm)
- endif
-
-! ***************************************************
-! end of Cuthill McKee permutation
-! ***************************************************
-
call get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declar.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declar.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declar.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -1,9 +1,4 @@
-!! DK DK added this temporarily
-!!!!!!!!!!!!! integer, dimension(NSPEC_CRUST_MANTLE) :: perm,invperm
-!! DK DK suppressed for now because useless when CUTHILL-MCKEE is off
- integer, dimension(1) :: perm,invperm
-
!! DK DK added this for merged version
!! DK DK stored in single precision for merged version, check if it precise enough (probably yes)
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_1D_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_1D_buffers.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_1D_buffers.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -54,7 +54,7 @@
yread1D_leftxi_lefteta, yread1D_rightxi_lefteta, yread1D_leftxi_righteta, yread1D_rightxi_righteta, &
zread1D_leftxi_lefteta, zread1D_rightxi_lefteta, zread1D_leftxi_righteta, zread1D_rightxi_righteta
- integer nspec,myrank,nglob_ori,nglob,ipoin1D,iregion
+ integer nspec,myrank,nglob_ori,nglob,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER
logical iMPIcut_xi(2,nspec)
@@ -78,29 +78,6 @@
! MPI 1D buffer element numbering
integer ispeccount,npoin1D,ix,iy,iz
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- integer, dimension(:), allocatable :: ibool_selected
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-
-! allocate arrays for message buffers with maximum size
-! define maximum size for message buffers
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(ibool_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(xstore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ystore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(zstore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ind(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ninseg(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(iglob(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(locval(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ifseg(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(iwork(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(work(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- endif
-
! write the MPI buffers for the left and right edges of the slice
! and the position of the points to check that the buffers are fine
@@ -144,44 +121,20 @@
npoin1D = npoin1D + 1
!! DK DK added this for merged version
if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
!! DK DK added this for merged
- ibool1D_leftxi_lefteta(npoin1D) = ibool(ix,iy,iz,ispec)
- xread1D_leftxi_lefteta(npoin1D) = xstore(ix,iy,iz,ispec)
- yread1D_leftxi_lefteta(npoin1D) = ystore(ix,iy,iz,ispec)
- zread1D_leftxi_lefteta(npoin1D) = zstore(ix,iy,iz,ispec)
- endif
+ ibool1D_leftxi_lefteta(npoin1D) = ibool(ix,iy,iz,ispec)
+ xread1D_leftxi_lefteta(npoin1D) = xstore(ix,iy,iz,ispec)
+ yread1D_leftxi_lefteta(npoin1D) = ystore(ix,iy,iz,ispec)
+ zread1D_leftxi_lefteta(npoin1D) = zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged version
- if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
-
- do ipoin1D=1,npoin1D
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
-!! DK DK suppressed merged ystore_selected(ipoin1D),zstore_selected(ipoin1D)
-!! DK DK added this for merged
- ibool1D_leftxi_lefteta(ipoin1D) = ibool_selected(ipoin1D)
- xread1D_leftxi_lefteta(ipoin1D) = xstore_selected(ipoin1D)
- yread1D_leftxi_lefteta(ipoin1D) = ystore_selected(ipoin1D)
- zread1D_leftxi_lefteta(ipoin1D) = zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -229,44 +182,20 @@
npoin1D = npoin1D + 1
!! DK DK added this for merged version
if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
!! DK DK added this for merged
- ibool1D_rightxi_lefteta(npoin1D) = ibool(ix,iy,iz,ispec)
- xread1D_rightxi_lefteta(npoin1D) = xstore(ix,iy,iz,ispec)
- yread1D_rightxi_lefteta(npoin1D) = ystore(ix,iy,iz,ispec)
- zread1D_rightxi_lefteta(npoin1D) = zstore(ix,iy,iz,ispec)
- endif
+ ibool1D_rightxi_lefteta(npoin1D) = ibool(ix,iy,iz,ispec)
+ xread1D_rightxi_lefteta(npoin1D) = xstore(ix,iy,iz,ispec)
+ yread1D_rightxi_lefteta(npoin1D) = ystore(ix,iy,iz,ispec)
+ zread1D_rightxi_lefteta(npoin1D) = zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged version
- if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
-
- do ipoin1D=1,npoin1D
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
-!! DK DK suppressed merged ystore_selected(ipoin1D),zstore_selected(ipoin1D)
-!! DK DK added this for merged
- ibool1D_rightxi_lefteta(ipoin1D) = ibool_selected(ipoin1D)
- xread1D_rightxi_lefteta(ipoin1D) = xstore_selected(ipoin1D)
- yread1D_rightxi_lefteta(ipoin1D) = ystore_selected(ipoin1D)
- zread1D_rightxi_lefteta(ipoin1D) = zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -324,44 +253,20 @@
npoin1D = npoin1D + 1
!! DK DK added this for merged version
if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
!! DK DK added this for merged
- ibool1D_leftxi_righteta(npoin1D) = ibool(ix,iy,iz,ispec)
- xread1D_leftxi_righteta(npoin1D) = xstore(ix,iy,iz,ispec)
- yread1D_leftxi_righteta(npoin1D) = ystore(ix,iy,iz,ispec)
- zread1D_leftxi_righteta(npoin1D) = zstore(ix,iy,iz,ispec)
- endif
+ ibool1D_leftxi_righteta(npoin1D) = ibool(ix,iy,iz,ispec)
+ xread1D_leftxi_righteta(npoin1D) = xstore(ix,iy,iz,ispec)
+ yread1D_leftxi_righteta(npoin1D) = ystore(ix,iy,iz,ispec)
+ zread1D_leftxi_righteta(npoin1D) = zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged version
- if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
-
- do ipoin1D=1,npoin1D
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
-!! DK DK suppressed merged ystore_selected(ipoin1D),zstore_selected(ipoin1D)
-!! DK DK added this for merged
- ibool1D_leftxi_righteta(ipoin1D) = ibool_selected(ipoin1D)
- xread1D_leftxi_righteta(ipoin1D) = xstore_selected(ipoin1D)
- yread1D_leftxi_righteta(ipoin1D) = ystore_selected(ipoin1D)
- zread1D_leftxi_righteta(ipoin1D) = zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -415,43 +320,19 @@
npoin1D = npoin1D + 1
!! DK DK added this for merged version
if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- ibool1D_rightxi_righteta(npoin1D) = ibool(ix,iy,iz,ispec)
- xread1D_rightxi_righteta(npoin1D) = xstore(ix,iy,iz,ispec)
- yread1D_rightxi_righteta(npoin1D) = ystore(ix,iy,iz,ispec)
- zread1D_rightxi_righteta(npoin1D) = zstore(ix,iy,iz,ispec)
- endif
+ ibool1D_rightxi_righteta(npoin1D) = ibool(ix,iy,iz,ispec)
+ xread1D_rightxi_righteta(npoin1D) = xstore(ix,iy,iz,ispec)
+ yread1D_rightxi_righteta(npoin1D) = ystore(ix,iy,iz,ispec)
+ zread1D_rightxi_righteta(npoin1D) = zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged version
- if(npoin1D > NGLOB1D_RADIAL_MAX) stop 'DK DK error merged NGLOB1D_RADIAL_MAX'
-
- do ipoin1D=1,npoin1D
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
-!! DK DK suppressed merged ystore_selected(ipoin1D),zstore_selected(ipoin1D)
-!! DK DK added this for merged
- ibool1D_rightxi_righteta(ipoin1D) = ibool_selected(ipoin1D)
- xread1D_rightxi_righteta(ipoin1D) = xstore_selected(ipoin1D)
- yread1D_rightxi_righteta(ipoin1D) = ystore_selected(ipoin1D)
- zread1D_rightxi_righteta(ipoin1D) = zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -464,19 +345,5 @@
if(ispeccount /= NSPEC1D_RADIAL_CORNER(iregion,3) .or. npoin1D /= NGLOB1D_RADIAL_CORNER(iregion,3)) &
call exit_MPI(myrank,'error MPI 1D buffer detection in xi=right')
- if (PERFORM_CUTHILL_MCKEE) then
- deallocate(ibool_selected)
- deallocate(xstore_selected)
- deallocate(ystore_selected)
- deallocate(zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
- endif
-
end subroutine get_MPI_1D_buffers
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_eta.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_eta.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -26,8 +26,8 @@
!=====================================================================
subroutine get_MPI_cutplanes_eta(myrank,nspec,iMPIcut_eta,ibool, &
- xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_XI_FACE,iregion,NGLOB2DMAX_XY,nglob_ori,iboolleft_eta,iboolright_eta,NGLOB2DMAX_YMIN_YMAX,npoin2D_eta)
+ mask_ibool,npointot, &
+ NSPEC2D_XI_FACE,iregion,nglob_ori,iboolleft_eta,iboolright_eta,NGLOB2DMAX_YMIN_YMAX,npoin2D_eta)
! this routine detects cut planes along eta
! In principle the left cut plane of the first slice
@@ -41,17 +41,13 @@
integer :: NGLOB2DMAX_YMIN_YMAX
integer, dimension(NGLOB2DMAX_YMIN_YMAX) :: iboolleft_eta,iboolright_eta
- integer nspec,myrank,nglob_ori,nglob,ipoin2D,NGLOB2DMAX_XY,iregion
+ integer nspec,myrank,nglob_ori,nglob,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_XI_FACE
logical iMPIcut_eta(2,nspec)
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
! logical mask used to create arrays iboolleft_eta and iboolright_eta
integer npointot
logical mask_ibool(npointot)
@@ -63,30 +59,6 @@
integer ispecc1,ispecc2,npoin2D_eta,ix,iy,iz
integer nspec2Dtheor
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- integer, dimension(:), allocatable :: ibool_selected
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-
-
-! allocate arrays for message buffers with maximum size
-! define maximum size for message buffers
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(ibool_selected(NGLOB2DMAX_XY))
- allocate(xstore_selected(NGLOB2DMAX_XY))
- allocate(ystore_selected(NGLOB2DMAX_XY))
- allocate(zstore_selected(NGLOB2DMAX_XY))
- allocate(ind(NGLOB2DMAX_XY))
- allocate(ninseg(NGLOB2DMAX_XY))
- allocate(iglob(NGLOB2DMAX_XY))
- allocate(locval(NGLOB2DMAX_XY))
- allocate(ifseg(NGLOB2DMAX_XY))
- allocate(iwork(NGLOB2DMAX_XY))
- allocate(work(NGLOB2DMAX_XY))
- endif
-
! theoretical number of surface elements in the buffers
! cut planes along eta=constant correspond to XI faces
nspec2Dtheor = NSPEC2D_XI_FACE(iregion,1)
@@ -123,19 +95,12 @@
npoin2D_eta = npoin2D_eta + 1
!! DK DK added this for merged
if(npoin2D_eta > NGLOB2DMAX_YMIN_YMAX) stop 'DK DK error points merged'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_eta) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_eta) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_eta) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_eta) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
!! DK DK added this for merged
!! DK DK merged ces deux tableaux sont les memes donc on pourrait n'en declarer qu'un seul
!! DK DK merged mais en fait non car on le reutilise ci-dessous pour ibool_right
- iboolleft_eta(npoin2D_eta) = ibool(ix,iy,iz,ispec)
- endif
+ iboolleft_eta(npoin2D_eta) = ibool(ix,iy,iz,ispec)
endif
enddo
enddo
@@ -143,23 +108,7 @@
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_eta,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged
- if(npoin2D_eta > NGLOB2DMAX_YMIN_YMAX) stop 'DK DK error points merged'
-
- do ipoin2D=1,npoin2D_eta
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
-!! DK DK suppressed merged ystore_selected(ipoin2D),zstore_selected(ipoin2D)
-!! DK DK added this for merged
-!! DK DK merged ces deux tableaux sont les memes donc on pourrait n'en declarer qu'un seul
-!! DK DK merged mais en fait non car on le reutilise ci-dessous pour ibool_right
- iboolleft_eta(ipoin2D) = ibool_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -201,19 +150,12 @@
npoin2D_eta = npoin2D_eta + 1
!! DK DK added this for merged
if(npoin2D_eta > NGLOB2DMAX_YMIN_YMAX) stop 'DK DK error points merged'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_eta) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_eta) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_eta) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_eta) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
!! DK DK added this for merged
!! DK DK merged ces deux tableaux sont les memes donc on pourrait n'en declarer qu'un seul
!! DK DK merged mais en fait non car on le reutilise ci-dessous pour ibool_right
- iboolright_eta(npoin2D_eta) = ibool(ix,iy,iz,ispec)
- endif
+ iboolright_eta(npoin2D_eta) = ibool(ix,iy,iz,ispec)
endif
enddo
enddo
@@ -221,23 +163,7 @@
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_eta,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged
- if(npoin2D_eta > NGLOB2DMAX_YMIN_YMAX) stop 'DK DK error points merged'
-
- do ipoin2D=1,npoin2D_eta
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
-!! DK DK suppressed merged ystore_selected(ipoin2D),zstore_selected(ipoin2D)
-!! DK DK added this for merged
-!! DK DK merged ces deux tableaux sont les memes donc on pourrait n'en declarer qu'un seul
-!! DK DK merged mais en fait non car on le reutilise ci-dessous pour ibool_right
- iboolright_eta(ipoin2D) = ibool_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -249,19 +175,5 @@
! compare number of surface elements detected to analytical value
if(ispecc2 /= nspec2Dtheor) call exit_MPI(myrank,'error MPI cut-planes detection in eta=right')
- if (PERFORM_CUTHILL_MCKEE) then
- deallocate(ibool_selected)
- deallocate(xstore_selected)
- deallocate(ystore_selected)
- deallocate(zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
- endif
-
end subroutine get_MPI_cutplanes_eta
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_xi.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_xi.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_MPI_cutplanes_xi.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -26,8 +26,8 @@
!=====================================================================
subroutine get_MPI_cutplanes_xi(myrank,nspec,iMPIcut_xi,ibool, &
- xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_ETA_FACE,iregion,NGLOB2DMAX_XY,nglob_ori,iboolleft_xi,iboolright_xi,NGLOB2DMAX_XMIN_XMAX,npoin2D_xi)
+ mask_ibool,npointot, &
+ NSPEC2D_ETA_FACE,iregion,nglob_ori,iboolleft_xi,iboolright_xi,NGLOB2DMAX_XMIN_XMAX,npoin2D_xi)
! this routine detects cut planes along xi
! In principle the left cut plane of the first slice
@@ -41,17 +41,13 @@
integer :: NGLOB2DMAX_XMIN_XMAX
integer, dimension(NGLOB2DMAX_XMIN_XMAX) :: iboolleft_xi,iboolright_xi
- integer nspec,myrank,nglob_ori,nglob,ipoin2D,iregion
+ integer nspec,myrank,nglob_ori,nglob,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_ETA_FACE
logical iMPIcut_xi(2,nspec)
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
! logical mask used to create arrays iboolleft_xi and iboolright_xi
integer npointot
logical mask_ibool(npointot)
@@ -65,31 +61,6 @@
character(len=150) errmsg
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- integer NGLOB2DMAX_XY
- integer, dimension(:), allocatable :: ibool_selected
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-
-! allocate arrays for message buffers with maximum size
-! define maximum size for message buffers
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(ibool_selected(NGLOB2DMAX_XY))
- allocate(xstore_selected(NGLOB2DMAX_XY))
- allocate(ystore_selected(NGLOB2DMAX_XY))
- allocate(zstore_selected(NGLOB2DMAX_XY))
- allocate(ind(NGLOB2DMAX_XY))
- allocate(ninseg(NGLOB2DMAX_XY))
- allocate(iglob(NGLOB2DMAX_XY))
- allocate(locval(NGLOB2DMAX_XY))
- allocate(ifseg(NGLOB2DMAX_XY))
- allocate(iwork(NGLOB2DMAX_XY))
- allocate(work(NGLOB2DMAX_XY))
- endif
-
-
! theoretical number of surface elements in the buffers
! cut planes along xi=constant correspond to ETA faces
nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,1)
@@ -125,17 +96,10 @@
npoin2D_xi = npoin2D_xi + 1
!! DK DK added this for merged
if(npoin2D_xi > NGLOB2DMAX_XMIN_XMAX) stop 'DK DK error points merged'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_xi) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_xi) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_xi) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_xi) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
!! DK DK added this for merged
- iboolleft_xi(npoin2D_xi) = ibool(ix,iy,iz,ispec)
- endif
+ iboolleft_xi(npoin2D_xi) = ibool(ix,iy,iz,ispec)
endif
enddo
enddo
@@ -143,23 +107,7 @@
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_xi,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged
- if(npoin2D_xi > NGLOB2DMAX_XMIN_XMAX) stop 'DK DK error points merged'
-
- do ipoin2D=1,npoin2D_xi
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
-!! DK DK suppressed merged ystore_selected(ipoin2D),zstore_selected(ipoin2D)
-!! DK DK added this for merged
-!! DK DK merged ces deux tableaux sont les memes donc on pourrait n'en declarer qu'un seul
-!! DK DK merged mais en fait non car on le reutilise ci-dessous pour ibool_right
- iboolleft_xi(ipoin2D) = ibool_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -203,16 +151,9 @@
npoin2D_xi = npoin2D_xi + 1
!! DK DK added this for merged
if(npoin2D_xi > NGLOB2DMAX_XMIN_XMAX) stop 'DK DK error points merged'
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_xi) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_xi) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_xi) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_xi) = zstore(ix,iy,iz,ispec)
- else
!! DK DK suppressed merged write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
!! DK DK suppressed merged ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- iboolright_xi(npoin2D_xi) = ibool(ix,iy,iz,ispec)
- endif
+ iboolright_xi(npoin2D_xi) = ibool(ix,iy,iz,ispec)
endif
enddo
enddo
@@ -220,20 +161,7 @@
enddo
nglob=nglob_ori
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_xi,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-!! DK DK added this for merged
- if(npoin2D_xi > NGLOB2DMAX_XMIN_XMAX) stop 'DK DK error points merged'
-
- do ipoin2D=1,npoin2D_xi
-!! DK DK suppressed merged write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
-!! DK DK suppressed merged ystore_selected(ipoin2D),zstore_selected(ipoin2D)
- iboolright_xi(ipoin2D) = ibool_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
!! DK DK suppressed merged write(10,*) '0 0 0. 0. 0.'
@@ -248,19 +176,5 @@
call exit_MPI(myrank,errmsg)
endif
- if (PERFORM_CUTHILL_MCKEE) then
- deallocate(ibool_selected)
- deallocate(xstore_selected)
- deallocate(ystore_selected)
- deallocate(zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
- endif
-
end subroutine get_MPI_cutplanes_xi
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_boundaries.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_boundaries.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -137,10 +137,6 @@
call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmin,normal_xmin,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb1,xstore,ystore,zstore,ibelm_xmin,normal_xmin,&
- jacobian2D_xmin,NSPEC2DMAX_XMIN_XMAX,NGLLY,NGLLZ,nspec)
- endif
endif
! on boundary: xmax
@@ -182,10 +178,6 @@
call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmax,normal_xmax,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb2,xstore,ystore,zstore,ibelm_xmax,normal_xmax,&
- jacobian2D_xmax,NSPEC2DMAX_XMIN_XMAX,NGLLY,NGLLZ,nspec)
- endif
endif
! on boundary: ymin
@@ -227,10 +219,6 @@
call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymin,normal_ymin,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb3,xstore,ystore,zstore,ibelm_ymin,normal_ymin,&
- jacobian2D_ymin,NSPEC2DMAX_YMIN_YMAX,NGLLX,NGLLZ,nspec)
- endif
endif
! on boundary: ymax
@@ -272,10 +260,6 @@
call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymax,normal_ymax,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb4,xstore,ystore,zstore,ibelm_ymax,normal_ymax,&
- jacobian2D_ymax,NSPEC2DMAX_YMIN_YMAX,NGLLX,NGLLZ,nspec)
- endif
endif
! on boundary: bottom
@@ -316,10 +300,6 @@
call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, &
jacobian2D_bottom,normal_bottom,NGLLX,NGLLY,NSPEC2D_BOTTOM)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb5,xstore,ystore,zstore,ibelm_bottom,normal_bottom,&
- jacobian2D_bottom,NSPEC2D_BOTTOM,NGLLX,NGLLY,nspec)
- endif
endif
! on boundary: top
@@ -360,10 +340,6 @@
call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb6,xstore,ystore,zstore,ibelm_top,normal_top,&
- jacobian2D_top,NSPEC2D_TOP,NGLLX,NGLLY,nspec)
- endif
endif
enddo
@@ -451,91 +427,3 @@
end subroutine compute_jacobian_2D
-
-
-subroutine sort_arrays_for_cuthill (ispecb,xstore,ystore,zstore,ibelm,normal,jacobian2D,nspec2D,NGLL1,NGLL2,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer :: ispecb,nspec2D,NGLL1,NGLL2,nspec,ispec_tmp,dummy_var,i
-
- integer ibelm(nspec2D)
- real(kind=CUSTOM_REAL) jacobian2D(NGLL1,NGLL2,NSPEC2D)
- real(kind=CUSTOM_REAL) normal(NDIM,NGLL1,NGLL2,NSPEC2D)
-
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
- integer, dimension(:), allocatable :: perm
- integer, dimension(:), allocatable :: ibelm_tmp
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_tmp
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_tmp
-
-! get permutation
- allocate (xstore_selected(ispecb))
- allocate (ystore_selected(ispecb))
- allocate (zstore_selected(ispecb))
- allocate(ind(ispecb))
- allocate(ninseg(ispecb))
- allocate(iglob(ispecb))
- allocate(locval(ispecb))
- allocate(ifseg(ispecb))
- allocate(iwork(ispecb))
- allocate(work(ispecb))
- allocate(perm(ispecb))
-
- do ispec_tmp=1,ispecb
- xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm(ispec_tmp))
- ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm(ispec_tmp))
- zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm(ispec_tmp))
- perm(ispec_tmp) = ispec_tmp
- enddo
-
- call sort_array_coordinates(ispecb,xstore_selected,ystore_selected,zstore_selected, &
- perm,iglob,locval,ifseg,dummy_var,ind,ninseg,iwork,work)
-
- deallocate (xstore_selected)
- deallocate (ystore_selected)
- deallocate (zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
-
-! permutation of ibelm
- allocate(ibelm_tmp(ispecb))
- ibelm_tmp(1:ispecb) = ibelm(1:ispecb)
- do i = 1,ispecb
- ibelm(perm(i)) = ibelm_tmp(i)
- enddo
- deallocate(ibelm_tmp)
-
-! permutation of normal
- allocate(normal_tmp(NDIM,NGLL1,NGLL2,ispecb))
- normal_tmp(:,:,:,1:ispecb) = normal(:,:,:,1:ispecb)
- do i = 1,ispecb
- normal(:,:,:,perm(i)) = normal_tmp(:,:,:,i)
- enddo
- deallocate(normal_tmp)
-
-! permutation of jacobian2D
- allocate(jacobian2D_tmp(NGLL1,NGLL2,ispecb))
- jacobian2D_tmp(:,:,1:ispecb) = jacobian2D(:,:,1:ispecb)
- do i = 1,ispecb
- jacobian2D(:,:,perm(i)) = jacobian2D_tmp(:,:,i)
- enddo
- deallocate(jacobian2D_tmp)
- deallocate(perm)
-
-end subroutine sort_arrays_for_cuthill
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_perm_cuthill_mckee.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_perm_cuthill_mckee.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_perm_cuthill_mckee.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -1,804 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory, California Institute of Technology, USA
-! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau / CNRS / INRIA
-! February 2008
-!
-! 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.
-!
-!=====================================================================
-
-! implement reverse Cuthill-McKee (1969) ordering, introduced in
-! E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices.
-! In Proceedings of the 1969 24th national conference, pages 157-172,
-! New-York, New-York, USA, 1969. ACM Press.
-! see for instance http://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm
-
- subroutine get_perm(ibool,perm,limit,nspec,nglob,INVERSE,FACE)
-
- implicit none
-
- include "constants.h"
-
- logical :: INVERSE,FACE
-
-! input
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! output
- integer, dimension(nspec) :: perm
-
-! local variables
- integer nspec,nglob_GLL_full
-
-! a neighbor of a hexahedral node is a hexahedra which share a face whith it -> max degre of a node = 6
- integer, parameter :: MAX_NUMBER_OF_NEIGHBORS = 100
-
-! global corner numbers that need to be created
- integer, dimension(nglob) :: global_corner_number
-
- integer mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
- integer, dimension(:), allocatable :: ne,np,adj
- integer xadj(nspec+1)
-
-! arrays to store the permutation and inverse permutation of the Cuthill-McKee algorithm
- integer, dimension(nspec) :: invperm
-
- logical maskel(nspec)
-
- integer i,istart,istop,number_of_neighbors
-
- integer nglob_eight_corners_only,nglob
-
-! only count the total size of the array that will be created, or actually create it
- logical count_only
- integer total_size_ne,total_size_adj,limit
-
-!
-!-----------------------------------------------------------------------
-!
- if(PERFORM_CUTHILL_MCKEE) then
-
- ! total number of points in the mesh
- nglob_GLL_full = nglob
-
- !---- call Charbel Farhat's routines
- call form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number,nglob_GLL_full,ibool,nglob_eight_corners_only)
- do i=1,nspec
- istart = mp(i)
- istop = mp(i+1) - 1
- enddo
-
- allocate(np(nglob_eight_corners_only+1))
- count_only = .true.
- total_size_ne = 1
- allocate(ne(total_size_ne))
- call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
- deallocate(ne)
- allocate(ne(total_size_ne))
- count_only = .false.
- call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
- do i=1,nglob_eight_corners_only
- istart = np(i)
- istop = np(i+1) - 1
- enddo
-
- count_only = .true.
- total_size_adj = 1
- allocate(adj(total_size_adj))
- call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
- count_only,total_size_ne,total_size_adj,FACE)
- deallocate(adj)
- allocate(adj(total_size_adj))
- count_only = .false.
- call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
- count_only,total_size_ne,total_size_adj,FACE)
- do i=1,nspec
- istart = xadj(i)
- istop = xadj(i+1) - 1
- number_of_neighbors = istop-istart+1
- if(number_of_neighbors < 1 .or. number_of_neighbors > MAX_NUMBER_OF_NEIGHBORS) stop 'incorrect number of neighbors'
- enddo
- deallocate(ne,np)
-
-! call the Cuthill-McKee sorting algorithm
- call cuthill_mckee(adj,xadj,perm,invperm,nspec,total_size_adj,limit,INVERSE)
- deallocate(adj)
- else
-! create identity permutation in order to do nothing
- do i=1,nspec
- perm(i) = i
- enddo
- endif
-
- end subroutine get_perm
-
-!=======================================================================
-!
-! Charbel Farhat's FEM topology routines
-!
-! Dimitri Komatitsch, February 1996 - Code based on Farhat's original version
-! described in his technical report from 1987
-!
-! modified and adapted by Dimitri Komatitsch, May 2006
-!
-!=======================================================================
-
- subroutine form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number, &
- nglob_GLL_full,ibool,nglob_eight_corners_only)
-
-!-----------------------------------------------------------------------
-!
-! Forms the MN and MP arrays
-!
-! Input :
-! -------
-! ibool Array needed to build the element connectivity table
-! nspec Number of elements in the domain
-! NGNOD_HEXAHEDRA number of nodes per hexahedron (brick with 8 corners)
-!
-! Output :
-! --------
-! MN, MP This is the element connectivity array pair.
-! Array MN contains the list of the element
-! connectivity, that is, the nodes contained in each
-! element. They are stored in a stacked fashion.
-!
-! Pointer array MP stores the location of each
-! element list. Its length is equal to the number
-! of elements plus one.
-!
-!-----------------------------------------------------------------------
-
- implicit none
-
- include "constants.h"
-
- integer nspec,nglob_GLL_full
-
-! arrays with mesh parameters per slice
- integer, intent(in), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! global corner numbers that need to be created
- integer, intent(out), dimension(nglob_GLL_full) :: global_corner_number
- integer, intent(out) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
- integer, intent(out) :: nglob_eight_corners_only
-
- integer ninter,nsum,ispec,node,k,inumcorner,ix,iy,iz
-
- ninter = 1
- nsum = 1
- mp(1) = 1
-
-!---- define topology of the elements in the mesh
-!---- we need to define adjacent numbers from the sub-mesh consisting of the corners only
- nglob_eight_corners_only = 0
- global_corner_number(:) = -1
-
- do ispec=1,nspec
-
- inumcorner = 0
- do iz = 1,NGLLZ,NGLLZ-1
- do iy = 1,NGLLY,NGLLY-1
- do ix = 1,NGLLX,NGLLX-1
-
- inumcorner = inumcorner + 1
- if(inumcorner > NGNOD_HEXAHEDRA) stop 'corner number too large'
-
-! check if this point was already assigned a number previously, otherwise create one and store it
- if(global_corner_number(ibool(ix,iy,iz,ispec)) == -1) then
- nglob_eight_corners_only = nglob_eight_corners_only + 1
- global_corner_number(ibool(ix,iy,iz,ispec)) = nglob_eight_corners_only
- endif
-
- node = global_corner_number(ibool(ix,iy,iz,ispec))
- do k=nsum,ninter-1
- if(node == mn(k)) goto 200
- enddo
-
- mn(ninter) = node
- ninter = ninter + 1
- 200 continue
-
- enddo
- enddo
- enddo
-
- nsum = ninter
- mp(ispec + 1) = nsum
-
- enddo
-
- end subroutine form_elt_connectivity_foelco
-
-!
-!----------------------------------------------------
-!
-
- subroutine form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only, &
- nspec,count_only,total_size_ne)
-
-!-----------------------------------------------------------------------
-!
-! Forms the NE and NP arrays
-!
-! Input :
-! -------
-! MN, MP, nspec
-! nglob_eight_corners_only Number of nodes in the domain
-!
-! Output :
-! --------
-! NE, NP This is the node-connected element array pair.
-! Integer array NE contains a list of the
-! elements connected to each node, stored in stacked fashion.
-!
-! Array NP is the pointer array for the
-! location of a node's element list in the NE array.
-! Its length is equal to the number of points plus one.
-!
-!-----------------------------------------------------------------------
-
- implicit none
-
- include "constants.h"
-
-! only count the total size of the array that will be created, or actually create it
- logical count_only
- integer total_size_ne
-
- integer nglob_eight_corners_only,nspec
-
- integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
-
- integer, intent(out) :: ne(total_size_ne),np(nglob_eight_corners_only+1)
-
- integer nsum,inode,ispec,j
-
- nsum = 1
- np(1) = 1
-
- do inode=1,nglob_eight_corners_only
- do 200 ispec=1,nspec
-
- do j=mp(ispec),mp(ispec + 1) - 1
- if (mn(j) == inode) then
- if(count_only) then
- total_size_ne = nsum
- else
- ne(nsum) = ispec
- endif
- nsum = nsum + 1
- goto 200
- endif
- enddo
- 200 continue
-
- np(inode + 1) = nsum
-
- enddo
-
- end subroutine form_node_connectivity_fonoco
-
-!
-!----------------------------------------------------
-!
-
- subroutine create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec, &
- nglob_eight_corners_only,count_only,total_size_ne,total_size_adj,face)
-
-!-----------------------------------------------------------------------
-!
-! Establishes the element adjacency information of the mesh
-! Two elements are considered adjacent if they share a face.
-!
-! Input :
-! -------
-! MN, MP, NE, NP, nspec
-! MASKEL logical mask (length = nspec)
-!
-! Output :
-! --------
-! ADJ, XADJ This is the element adjacency array pair. Array
-! ADJ contains the list of the elements adjacent to
-! element i. They are stored in a stacked fashion.
-! Pointer array XADJ stores the location of each element list.
-!
-!-----------------------------------------------------------------------
-
- implicit none
-
- include "constants.h"
-
-! only count the total size of the array that will be created, or actually create it
- logical count_only,face
- integer total_size_ne,total_size_adj
-
- integer nglob_eight_corners_only
-
- integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1),ne(total_size_ne),np(nglob_eight_corners_only+1)
-
- integer, intent(out) :: adj(total_size_adj),xadj(nspec+1)
-
- logical maskel(nspec)
- integer countel(nspec)
-
- integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel
-
- xadj(1) = 1
- iad = 1
-
- do ispec=1,nspec
-
-! reset mask
- maskel(:) = .false.
-
-! mask current element
- maskel(ispec) = .true.
- if (face) countel(:) = 0
-
- istart = mp(ispec)
- istop = mp(ispec+1) - 1
- do ino=istart,istop
- node = mn(ino)
- jstart = np(node)
- jstop = np(node + 1) - 1
- do 120 jel=jstart,jstop
- nelem = ne(jel)
- if(maskel(nelem)) goto 120
- if (face) then
- ! if 2 elements share at least 3 corners, therefore they share a face
- countel(nelem) = countel(nelem) + 1
- if (countel(nelem)>=3) then
- if(count_only) then
- total_size_adj = iad
- else
- adj(iad) = nelem
- endif
- maskel(nelem) = .true.
- iad = iad + 1
- endif
- else
- if(count_only) then
- total_size_adj = iad
- else
- adj(iad) = nelem
- endif
- maskel(nelem) = .true.
- iad = iad + 1
- endif
- 120 continue
- enddo
-
- xadj(ispec+1) = iad
-
- enddo
-
- end subroutine create_adjacency_table_adjncy
-
-!
-!----------------------------------------------------
-!
-
- subroutine cuthill_mckee(adj,xadj,mask,invperm_all,nspec,total_size_adj,limit,INVERSE)
-
- implicit none
- include "constants.h"
-
- integer, intent(in) :: nspec,total_size_adj, limit
- integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
-
- integer, intent(out), dimension(nspec) :: mask,invperm_all
- integer, dimension(nspec) :: invperm_sub
- logical :: INVERSE
- integer ispec,gsize,counter,nspec_sub,root,total_ordered_elts, next_root
-
-! fill the mask with ones
- mask(:) = 1
- invperm_all(:) = 0
- counter = 0
- nspec_sub = limit
- root = 1
- total_ordered_elts = 0
-
- do while(total_ordered_elts < nspec)
- ! creation of a sublist of sorted elements which fit in the cache (the criterion of size is limit)
- ! limit = nb of element that can fit in the L2 cache
- call Cut_McK( root, nspec, total_size_adj, xadj, adj, mask, gsize, invperm_sub, limit, nspec_sub, next_root)
- ! add the sublist in the main permutation list
- invperm_all(total_ordered_elts+1:total_ordered_elts+nspec_sub) = invperm_sub(1:nspec_sub)
- total_ordered_elts = total_ordered_elts + nspec_sub
- ! seek for a new root to build the new sublist
- if (next_root > 0) then
- root = next_root
- else
- if (total_ordered_elts /= nspec) &
- call find_next_root(next_root,xadj,adj,total_size_adj,mask,invperm_all,total_ordered_elts,nspec)
- root = next_root
- endif
- enddo
-
- if (INVERSE) then
- do ispec=1,nspec
- mask(invperm_all(ispec)) = ispec
- enddo
- else
- mask(:) = invperm_all(:)
- endif
-
- end subroutine cuthill_mckee
-
-
-!*******************************************************************************
-! Objective: Cuthill-McKee ordering
-! The algorithm is:
-!
-! X(1) = ROOT.
-! for ( I = 1 to N-1)
-! Find all unlabeled neighbors of X(I),
-! assign them the next available labels, in order of increasing degree.
-!
-! Parameters:
-! root the starting point for the cm ordering.
-! nbnodes the number of nodes.
-! nnz the number of adjacency entries.
-!
-! xadj/adj the graph
-! mask only those nodes with nonzero mask are considered
-!
-! gsize the number of the connected component
-! invp Inverse invputation (from new order to old order)
-!*******************************************************************************
-
-subroutine find_next_root(next_root,xadj,adj,total_size_adj,mask,invperm_all,total_ordered_elts,nspec)
-
- implicit none
-
- include "constants.h"
-
-! input
- integer, intent(in) :: total_size_adj,total_ordered_elts,nspec
- integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
- integer, intent(in), dimension(nspec) :: mask,invperm_all
-! output
- integer, intent(out) :: next_root
-! variables
- integer :: cur_node,neighbor_node,i,j
-
- do i=total_ordered_elts, 1, -1
- cur_node = invperm_all(i)
- do j= xadj(cur_node), xadj(cur_node+1)-1
- neighbor_node = adj(j)
- if (mask(neighbor_node)/=0) then
- next_root=neighbor_node
- return
- endif
- enddo
- enddo
-
-end subroutine find_next_root
-
-!*******************************************************************************
-! Objective: Cuthill-McKee ordering
-! The algorithm is:
-!
-! X(1) = ROOT.
-! for ( I = 1 to N-1)
-! Find all unlabeled neighbors of X(I),
-! assign them the next available labels, in order of increasing degree.
-!
-! Parameters:
-! root the starting point for the cm ordering.
-! nbnodes the number of nodes.
-! nnz the number of adjacency entries.
-!
-! xadj/adj the graph
-! mask only those nodes with nonzero mask are considered
-!
-! gsize the number of the connected component
-! invp Inverse invputation (from new order to old order)
-!*******************************************************************************
-
-subroutine Cut_McK( root, nbnodes, nnz, xadj, adj, mask, gsize, invp, limit, nspec_sub, next_root)
-
- implicit none
-
- include "constants.h"
-
-!--------------------------------------------------------------- Input Variables
- integer root, nnz, nbnodes, limit, nspec_sub, next_root
-
- integer xadj(nbnodes+1), adj(nnz), mask(nbnodes)
-
-!-------------------------------------------------------------- Output Variables
- integer gsize
- integer invp(nbnodes)
-
-!--------------------------------------------------------------- Local Variables
- integer i, j, k, l, lbegin, lnbr, linvp, lvlend, nbr, node, fnbr
- integer deg(nbnodes)
-
-! Find the degrees of the nodes in the subgraph specified by mask and root
-! Here invp is used to store a levelization of the subgraph
- invp(:)=0
- deg(:)=0
- call degree ( root, nbnodes, nnz, xadj, adj, mask, gsize, deg, invp)
-
- mask(root) = 0
-
- IF (gsize > 1) THEN
- !If there is at least 2 nodes in the subgraph
- lvlend = 0
- lnbr = 1
-
- DO while (lvlend < lnbr)
- !lbegin/lvlend point to the begin/end of the present level
- lbegin = lvlend + 1
- lvlend = lnbr
-
- do i= lbegin, lvlend
- node = invp(i)
-
- !Find the unnumbered neighbours of node.
- !fnbr/lnbr point to the first/last neighbors of node
- fnbr = lnbr + 1
- do j= xadj(node), xadj(node+1)-1
- nbr = adj(j)
-
- if (mask(nbr) /= 0) then
- lnbr = lnbr + 1
- mask(nbr) = 0
- invp(lnbr) = nbr
- endif
- enddo
-
- !If no neighbors, go to next node in this level.
- IF (lnbr > fnbr) THEN
- !Sort the neighbors of NODE in increasing order by degree.
- !Linear insertion is used.
- k = fnbr
- do while (k < lnbr)
- l = k
- k = k + 1
- nbr = invp(k)
-
- DO WHILE (fnbr < l)
- linvp = invp(l)
-
- if (deg(linvp) <= deg(nbr)) then
- exit
- endif
-
- invp(l+1) = linvp
- l = l-1
- ENDDO
-
- invp(l+1) = nbr
- enddo
- ENDIF
- enddo
- ENDDO
-
- ENDIF
-
- if (gsize > limit) then
- do i = limit + 1 , nbnodes
- node=invp(i)
- if (node /=0) mask(node) = 1
- enddo
- next_root = invp(limit +1)
- nspec_sub = limit
- else
- next_root = -1
- nspec_sub = gsize
- endif
-
-END subroutine Cut_McK
-
-
-!*******************************************************************************
-! Objective: computes the degrees of the nodes in the connected graph
-!
-! Parameters:
-! root the root node
-! nbnodes the number of nodes in the graph
-! nnz the graph size
-! xadj/adj the whole graph
-! mask Only nodes with mask == 0 are considered
-!
-! gsize the number of nodes in the connected graph
-! deg degree for all the nodes in the connected graph
-! level levelization of the connected graph
-!
-!*******************************************************************************
-
-subroutine degree( root, nbnodes, nnz, xadj, adj, mask, gsize, deg, level )
-
- implicit none
-
-!--------------------------------------------------------------- Input Variables
- integer root, nbnodes, nnz
- integer xadj(nbnodes+1), adj(nnz), mask(nbnodes)
-
-!-------------------------------------------------------------- Output Variables
- integer gsize
- integer deg(nbnodes), level(nbnodes)
-
-!--------------------------------------------------------------- Local Variables
- integer i, j, ideg, lbegin, lvlend, lvsize, nxt, nbr, node
-
-! The sign of xadj(I) is used to indicate if node i has been considered
- xadj(root) = -xadj(root)
- level(1) = root
- nxt = 1
- lvlend = 0
- lvsize = 1
-
- DO WHILE (lvsize > 0)
- ! lbegin/lvlend points the begin/end of the present level
- lbegin = lvlend + 1
- lvlend = nxt
-
- ! Find the degrees of nodes in the present level and generate the next level
- DO i= lbegin, lvlend
- node = level(i)
- ideg = 0
- do j= ABS( xadj(node) ), ABS( xadj(node+1) )-1
- nbr = adj(j)
-
- if (mask(nbr) /= 0) then
- ideg = ideg + 1
-
- if (xadj(nbr) >= 0) then
- xadj(nbr) = -xadj(nbr)
- nxt = nxt + 1
- level(nxt) = nbr
- endif
- endif
- enddo
-
- deg(node) = ideg
- ENDDO
-
- !Compute the level size of the next level
- lvsize = nxt - lvlend
- ENDDO
-
- !Reset xadj to its correct sign
- do i = 1, nxt
- node = level(i)
- xadj(node) = -xadj(node)
- enddo
-
- gsize = nxt
-
-END subroutine degree
-
-!
-!-----------------------------------------------------------------------
-!
-
- subroutine permute_elements_real(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- real(kind=CUSTOM_REAL), intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
- enddo
-
- end subroutine permute_elements_real
-
-!
-!-----------------------------------------------------------------------
-!
-
-!! DK DK added this for merged version
- subroutine permute_elements_xelm_yelm_zelm(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- real(kind=CUSTOM_REAL), intent(inout), dimension(NGNOD,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:) = array_to_permute(:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,new_ispec) = temp_array(:,old_ispec)
- enddo
-
- end subroutine permute_elements_xelm_yelm_zelm
-
-!
-!-----------------------------------------------------------------------
-!
-
-! implement permutation of elements for arrays of integer type
- subroutine permute_elements_integer(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- integer, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
- enddo
-
- end subroutine permute_elements_integer
-
-!
-!-----------------------------------------------------------------------
-!
-
-! implement permutation of elements for arrays of double precision type
- subroutine permute_elements_dble(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- double precision, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
- enddo
-
- end subroutine permute_elements_dble
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90 2008-08-06 14:12:52 UTC (rev 12544)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90 2008-08-06 14:41:41 UTC (rev 12545)
@@ -1566,7 +1566,6 @@
NSPEC2D_BOTTOM(iregion_code),NSPEC2D_TOP(iregion_code),ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE,NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
- max(NGLOB2DMAX_XMIN_XMAX(iregion_code),NGLOB2DMAX_YMIN_YMAX(iregion_code)), &
myrank,LOCAL_PATH,OCEANS,ibathy_topo,rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD, &
ATTENUATION,ATTENUATION_3D,NCHUNKS,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
@@ -1593,7 +1592,7 @@
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle, &
rmass_crust_mantle,xelm_store_crust_mantle,yelm_store_crust_mantle,zelm_store_crust_mantle, &
!! DK DK this will have to change to fully support David's code to cut the superbrick
- npoin2D_xi_crust_mantle(1),npoin2D_eta_crust_mantle(1),perm,invperm)
+ npoin2D_xi_crust_mantle(1),npoin2D_eta_crust_mantle(1))
else if(iregion_code == IREGION_OUTER_CORE) then
! outer_core
@@ -1605,7 +1604,6 @@
ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE,NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
- max(NGLOB2DMAX_XMIN_XMAX(iregion_code),NGLOB2DMAX_YMIN_YMAX(iregion_code)), &
myrank,LOCAL_PATH,OCEANS,ibathy_topo,rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD, &
ATTENUATION,ATTENUATION_3D,NCHUNKS,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
@@ -1632,7 +1630,7 @@
kappavstore_outer_core,kappahstore_outer_core,muvstore_outer_core,muhstore_outer_core,eta_anisostore_outer_core, &
rmass_outer_core,xelm_store_outer_core,yelm_store_outer_core,zelm_store_outer_core, &
!! DK DK this will have to change to fully support David's code to cut the superbrick
- npoin2D_xi_outer_core(1),npoin2D_eta_outer_core(1),perm,invperm)
+ npoin2D_xi_outer_core(1),npoin2D_eta_outer_core(1))
else if(iregion_code == IREGION_INNER_CORE) then
! inner_core
@@ -1644,7 +1642,6 @@
ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE,NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
- max(NGLOB2DMAX_XMIN_XMAX(iregion_code),NGLOB2DMAX_YMIN_YMAX(iregion_code)), &
myrank,LOCAL_PATH,OCEANS,ibathy_topo,rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD, &
ATTENUATION,ATTENUATION_3D,NCHUNKS,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
@@ -1670,7 +1667,7 @@
kappavstore_inner_core,kappahstore_inner_core,muvstore_inner_core,muhstore_inner_core,eta_anisostore_inner_core, &
rmass_inner_core,xelm_store_inner_core,yelm_store_inner_core,zelm_store_inner_core, &
!! DK DK this will have to change to fully support David's code to cut the superbrick
- npoin2D_xi_inner_core(1),npoin2D_eta_inner_core(1),perm,invperm)
+ npoin2D_xi_inner_core(1),npoin2D_eta_inner_core(1))
else
stop 'DK DK incorrect region in merged code'
More information about the cig-commits
mailing list