From 5db0729e61ae3ec986aef064aa183a670da81eb7 Mon Sep 17 00:00:00 2001 From: Dimitri Komatitsch Date: Wed, 19 Feb 2014 01:34:57 +0100 Subject: [PATCH] suppressed ipass, PERFORM_CUTHILL_MCKEE, ACTUALLY_IMPLEMENT_PERM* and FURTHER_REDUCE_CACHE_MISSES, which were making the code too complex with no significant performance improvement. Also suppressed GENERATE_PARAVER_TRACES, which is not used any more. --- setup/constants.h.in | 21 - src/meshfem2D/read_parameter_file.F90 | 6 - src/meshfem2D/save_databases.f90 | 3 - src/specfem2D/Makefile.in | 4 - src/specfem2D/createnum_fast.f90 | 12 +- src/specfem2D/createnum_slow.f90 | 15 +- src/specfem2D/get_MPI.F90 | 12 +- src/specfem2D/get_global.f90 | 47 +- src/specfem2D/get_perm_cuthill_mckee.f90 | 726 ------------------ src/specfem2D/gmat01.f90 | 8 +- src/specfem2D/initialize_simulation.F90 | 14 +- src/specfem2D/locate_receivers.F90 | 8 +- src/specfem2D/locate_source_force.F90 | 8 +- src/specfem2D/locate_source_moment_tensor.F90 | 11 +- src/specfem2D/read_databases.F90 | 137 +--- src/specfem2D/set_sources.f90 | 23 +- src/specfem2D/setup_sources_receivers.F90 | 20 +- src/specfem2D/specfem2D.F90 | 330 ++------ 18 files changed, 147 insertions(+), 1258 deletions(-) delete mode 100644 src/specfem2D/get_perm_cuthill_mckee.f90 diff --git a/setup/constants.h.in b/setup/constants.h.in index b03d6b543..4f945c711 100644 --- a/setup/constants.h.in +++ b/setup/constants.h.in @@ -33,33 +33,12 @@ ! maximum size of the X and Z directions of a JPEG image generated by the display routine integer, parameter :: NX_NZ_IMAGE_MAX = 10000 -! further reduce cache misses inner/outer in two passes in the case of an MPI simulation -! this flag is ignored in the case of a serial simulation - logical, parameter :: FURTHER_REDUCE_CACHE_MISSES = .true. - -! for inverse Cuthill-McKee (1969) permutation - logical, parameter :: INVERSE = .true. - logical, parameter :: FACE = .false. - integer, parameter :: NGNOD_QUADRANGLE = 4 -! 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 - -! implement Cuthill-McKee or replace with identity permutation - logical, parameter :: ACTUALLY_IMPLEMENT_PERM_OUT = .false. - logical, parameter :: ACTUALLY_IMPLEMENT_PERM_INN = .false. - logical, parameter :: ACTUALLY_IMPLEMENT_PERM_WHOLE = .true. - ! create file DATA/model_velocity.dat_output or not (can be huge and slow, thus off by default) logical, parameter :: OUTPUT_MODEL_VELOCITY_FILE = .false. ! in PostScript files about mesh quality, show all elements above this threshold double precision, parameter :: THRESHOLD_POSTSCRIPT = 95.d0 / 100.d0 -! add MPI barriers and suppress seismograms if we generate traces of the run for analysis with "ParaVer" - logical, parameter :: GENERATE_PARAVER_TRACES = .false. - ! option to display only part of the mesh and not the whole mesh, ! for instance to analyze Cuthill-McKee mesh partitioning etc. ! Possible values are: diff --git a/src/meshfem2D/read_parameter_file.F90 b/src/meshfem2D/read_parameter_file.F90 index b1ae89f1e..7ca3fe046 100644 --- a/src/meshfem2D/read_parameter_file.F90 +++ b/src/meshfem2D/read_parameter_file.F90 @@ -137,9 +137,6 @@ module parameter_file ! non linear display to enhance small amplitudes in color images double precision :: POWER_DISPLAY_COLOR -! perform inverse Cuthill-McKee (1969) permutation for mesh numbering - logical :: PERFORM_CUTHILL_MCKEE - ! output seismograms in Seismic Unix format (adjoint traces will be read in the same format) logical :: SU_FORMAT @@ -201,9 +198,6 @@ subroutine read_parameter_file() call read_value_integer_p(partitioning_method, 'mesher.partitioning_method') if(err_occurred() /= 0) stop 'error reading parameter 5a in Par_file' - call read_value_logical_p(PERFORM_CUTHILL_MCKEE, 'mesher.PERFORM_CUTHILL_MCKEE') - if(err_occurred() /= 0) stop 'error reading parameter 5b in Par_file' - call read_value_integer_p(ngnod, 'mesher.ngnod') if(err_occurred() /= 0) stop 'error reading parameter 6 in Par_file' diff --git a/src/meshfem2D/save_databases.f90 b/src/meshfem2D/save_databases.f90 index e2084b78a..c09bfb20b 100644 --- a/src/meshfem2D/save_databases.f90 +++ b/src/meshfem2D/save_databases.f90 @@ -197,9 +197,6 @@ subroutine save_databases(nspec,num_material,region_pml_external_mesh, & write(15,*) 'POWER_DISPLAY_COLOR' write(15,*) POWER_DISPLAY_COLOR - write(15,*) 'PERFORM_CUTHILL_MCKEE' - write(15,*) PERFORM_CUTHILL_MCKEE - write(15,*) 'SU_FORMAT' write(15,*) SU_FORMAT diff --git a/src/specfem2D/Makefile.in b/src/specfem2D/Makefile.in index e4bdf53cc..e68f856c4 100644 --- a/src/specfem2D/Makefile.in +++ b/src/specfem2D/Makefile.in @@ -128,7 +128,6 @@ OBJS_SPECFEM2D = \ $O/force_ftz.o \ $O/get_global.o \ $O/get_MPI.o \ - $O/get_perm_cuthill_mckee.o \ $O/get_poroelastic_velocities.o \ $O/gll_library.o \ $O/gmat01.o \ @@ -309,9 +308,6 @@ $O/get_global.o: ${S}/get_global.f90 ${SETUP}/constants.h $O/get_MPI.o: ${S}/get_MPI.F90 ${SETUP}/constants.h ${F90} $(FLAGS_CHECK) -c -o $O/get_MPI.o ${S}/get_MPI.F90 -$O/get_perm_cuthill_mckee.o: ${S}/get_perm_cuthill_mckee.f90 ${SETUP}/constants.h - ${F90} $(FLAGS_CHECK) -c -o $O/get_perm_cuthill_mckee.o ${S}/get_perm_cuthill_mckee.f90 - $O/get_poroelastic_velocities.o: ${S}/get_poroelastic_velocities.f90 ${SETUP}/constants.h ${F90} $(FLAGS_CHECK) -c -o $O/get_poroelastic_velocities.o ${S}/get_poroelastic_velocities.f90 diff --git a/src/specfem2D/createnum_fast.f90 b/src/specfem2D/createnum_fast.f90 index 3df81e0ec..56bd5e7ff 100644 --- a/src/specfem2D/createnum_fast.f90 +++ b/src/specfem2D/createnum_fast.f90 @@ -42,7 +42,7 @@ ! !======================================================================== - subroutine createnum_fast(knods,ibool,shape,coorg,nglob,npgeo,nspec,ngnod,myrank,ipass) + subroutine createnum_fast(knods,ibool,shape,coorg,nglob,npgeo,nspec,ngnod,myrank) ! same as subroutine "createnum_slow" but with a faster algorithm @@ -50,7 +50,7 @@ subroutine createnum_fast(knods,ibool,shape,coorg,nglob,npgeo,nspec,ngnod,myrank include "constants.h" - integer nglob,npgeo,nspec,ngnod,myrank,ipass + integer nglob,npgeo,nspec,ngnod,myrank integer knods(ngnod,nspec),ibool(NGLLX,NGLLZ,nspec) double precision shape(ngnod,NGLLX,NGLLX) double precision coorg(NDIM,npgeo) @@ -70,7 +70,7 @@ subroutine createnum_fast(knods,ibool,shape,coorg,nglob,npgeo,nspec,ngnod,myrank !---- create global mesh numbering - if(myrank == 0 .and. ipass == 1) then + if(myrank == 0) then write(IOUT,*) write(IOUT,*) write(IOUT,*) 'Generating global mesh numbering (fast version)...' @@ -231,12 +231,6 @@ subroutine createnum_fast(knods,ibool,shape,coorg,nglob,npgeo,nspec,ngnod,myrank ! check the numbering obtained if(minval(ibool) /= 1 .or. maxval(ibool) /= nglob) call exit_MPI('Error while generating global numbering') -! if(myrank == 0 .and. ipass == 1) then -! write(IOUT,*) -! write(IOUT,*) 'Total number of points of the global mesh on slice 0: ',nglob -! write(IOUT,*) -! endif - end subroutine createnum_fast diff --git a/src/specfem2D/createnum_slow.f90 b/src/specfem2D/createnum_slow.f90 index 74356940c..fd71351b9 100644 --- a/src/specfem2D/createnum_slow.f90 +++ b/src/specfem2D/createnum_slow.f90 @@ -42,7 +42,7 @@ ! !======================================================================== - subroutine createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank,ipass) + subroutine createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank) ! generate the global numbering @@ -50,7 +50,7 @@ subroutine createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank,ipass) include "constants.h" - integer nglob,nspec,ngnod,myrank,ipass + integer nglob,nspec,ngnod,myrank integer knods(ngnod,nspec),ibool(NGLLX,NGLLZ,nspec) @@ -63,7 +63,7 @@ subroutine createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank,ipass) !---- create global mesh numbering - if(myrank == 0 .and. ipass == 1) then + if(myrank == 0) then write(IOUT,*) write(IOUT,*) 'Generating global mesh numbering (slow version)...' write(IOUT,*) @@ -311,14 +311,5 @@ subroutine createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank,ipass) ! verification de la coherence de la numerotation generee if(minval(ibool) /= 1 .or. maxval(ibool) /= nglob) call exit_MPI('Error while generating global numbering') -! if(myrank == 0 .and. ipass == 1) then -! write(IOUT,*) 'Total number of points of the global mesh on slice 0: ',nglob,' distributed as follows:' -! write(IOUT,*) -! write(IOUT,*) 'Number of interior points: ',nglob-npedge-npcorn -! write(IOUT,*) 'Number of edge points (without corners): ',npedge -! write(IOUT,*) 'Number of corner points: ',npcorn -! write(IOUT,*) -! endif - end subroutine createnum_slow diff --git a/src/specfem2D/get_MPI.F90 b/src/specfem2D/get_MPI.F90 index 38e5d0fbe..63ca1def6 100644 --- a/src/specfem2D/get_MPI.F90 +++ b/src/specfem2D/get_MPI.F90 @@ -55,7 +55,7 @@ subroutine get_MPI(nspec,ibool,knods,ngnod,nglob,elastic,poroelastic, & inum_interfaces_poroelastic, & ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, & mask_ispec_inner_outer, & - myrank,ipass,coord) + myrank,coord) ! sets up the MPI interface for communication between partitions @@ -85,7 +85,7 @@ subroutine get_MPI(nspec,ibool,knods,ngnod,nglob,elastic,poroelastic, & logical, dimension(nspec), intent(inout) :: mask_ispec_inner_outer - integer :: myrank,ipass + integer :: myrank double precision, dimension(NDIM,nglob) :: coord !local parameters @@ -231,7 +231,7 @@ subroutine get_MPI(nspec,ibool,knods,ngnod,nglob,elastic,poroelastic, & ! outputs total number of MPI interface points call MPI_REDUCE(num_points2, num_points1, 1, MPI_INTEGER, & MPI_SUM, 0, MPI_COMM_WORLD, ier) - if( myrank == 0 .and. ipass == 1 ) then + if( myrank == 0 ) then write(IOUT,*) 'total MPI interface points: ',num_points1 endif @@ -276,7 +276,7 @@ subroutine get_MPI(nspec,ibool,knods,ngnod,nglob,elastic,poroelastic, & call MPI_REDUCE(inum, num_points1, 1, MPI_INTEGER, & MPI_SUM, 0, MPI_COMM_WORLD, ier) - if( myrank == 0 .and. ipass == 1 ) then + if( myrank == 0 ) then write(IOUT,*) ' acoustic interface points: ',num_points1 endif @@ -310,9 +310,7 @@ subroutine get_MPI(nspec,ibool,knods,ngnod,nglob,elastic,poroelastic, & MPI_SUM, 0, MPI_COMM_WORLD, ier) if( myrank == 0 ) then - if( ipass == 1 ) then - write(IOUT,*) ' assembly acoustic MPI interface points:',num_points2 - endif + write(IOUT,*) ' assembly acoustic MPI interface points:',num_points2 ! they don't need to fit, somehow.. !if( num_points2 /= num_points1 ) then diff --git a/src/specfem2D/get_global.f90 b/src/specfem2D/get_global.f90 index 8ce5332c2..a89840282 100644 --- a/src/specfem2D/get_global.f90 +++ b/src/specfem2D/get_global.f90 @@ -42,13 +42,12 @@ ! !======================================================================== - - subroutine get_global(nspec_outer,nspec,nglob,ibool) + subroutine get_global(nspec,nglob,ibool) implicit none include "constants.h" - integer :: nspec_outer,nspec,nglob + integer :: nspec,nglob integer, dimension(NGLLX,NGLLZ,nspec) :: ibool @@ -69,46 +68,6 @@ subroutine get_global(nspec_outer,nspec,nglob,ibool) inumber = 0 - if(.not. ACTUALLY_IMPLEMENT_PERM_WHOLE) then - - ! first reduce cache misses in outer elements, since they are taken first - ! loop over spectral elements - do ispec = 1,nspec_outer - do j=1,NGLLZ - do i=1,NGLLX - if(mask_ibool(copy_ibool_ori(i,j,ispec)) == -1) then - ! create a new point - inumber = inumber + 1 - ibool(i,j,ispec) = inumber - mask_ibool(copy_ibool_ori(i,j,ispec)) = inumber - else - ! use an existing point created previously - ibool(i,j,ispec) = mask_ibool(copy_ibool_ori(i,j,ispec)) - endif - enddo - enddo - enddo - - ! then reduce cache misses in inner elements, since they are taken second - ! loop over spectral elements - do ispec = nspec_outer+1,nspec - do j=1,NGLLZ - do i=1,NGLLX - if(mask_ibool(copy_ibool_ori(i,j,ispec)) == -1) then - ! create a new point - inumber = inumber + 1 - ibool(i,j,ispec) = inumber - mask_ibool(copy_ibool_ori(i,j,ispec)) = inumber - else - ! use an existing point created previously - ibool(i,j,ispec) = mask_ibool(copy_ibool_ori(i,j,ispec)) - endif - enddo - enddo - enddo - - else ! if ACTUALLY_IMPLEMENT_PERM_WHOLE - ! reduce cache misses in all the elements ! loop over spectral elements do ispec = 1,nspec @@ -127,8 +86,6 @@ subroutine get_global(nspec_outer,nspec,nglob,ibool) enddo enddo - endif - deallocate(mask_ibool,copy_ibool_ori) end subroutine get_global diff --git a/src/specfem2D/get_perm_cuthill_mckee.f90 b/src/specfem2D/get_perm_cuthill_mckee.f90 deleted file mode 100644 index 79506f9df..000000000 --- a/src/specfem2D/get_perm_cuthill_mckee.f90 +++ /dev/null @@ -1,726 +0,0 @@ - -!======================================================================== -! -! S P E C F E M 2 D Version 7 . 0 -! -------------------------------- -! -! Copyright CNRS, Inria and University of Pau, France, -! and Princeton University / California Institute of Technology, USA. -! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr -! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr -! Roland Martin, roland DOT martin aT univ-pau DOT fr -! Christina Morency, cmorency aT princeton DOT edu -! -! This software is a computer program whose purpose is to solve -! the two-dimensional viscoelastic anisotropic or poroelastic wave equation -! using a spectral-element method (SEM). -! -! This software is governed by the CeCILL license under French law and -! abiding by the rules of distribution of free software. You can use, -! modify and/or redistribute the software under the terms of the CeCILL -! license as circulated by CEA, CNRS and Inria at the following URL -! "http://www.cecill.info". -! -! As a counterpart to the access to the source code and rights to copy, -! modify and redistribute granted by the license, users are provided only -! with a limited warranty and the software's author, the holder of the -! economic rights, and the successive licensors have only limited -! liability. -! -! In this respect, the user's attention is drawn to the risks associated -! with loading, using, modifying and/or developing or reproducing the -! software by the user in light of its specific status of free software, -! that may mean that it is complicated to manipulate, and that also -! therefore means that it is reserved for developers and experienced -! professionals having in-depth computer knowledge. Users are therefore -! encouraged to load and test the software's suitability as regards their -! requirements in conditions enabling the security of their systems and/or -! data to be ensured and, more generally, to use and operate it in the -! same conditions as regards security. -! -! The full text of the license is available in file "LICENSE". -! -!======================================================================== - -! 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_cuthill_mckee(ibool,perm,limit,nspec,nglob,PERFORM_CUTHILL_MCKEE) - - implicit none - - include "constants.h" - -! local variables - integer nspec,nglob_GLL_full - integer nglob_four_corners_only,nglob - -! maximum number of neighbors of a spectral element (in principle, it could be any value) - integer, parameter :: MAX_NUMBER_OF_NEIGHBORS = 50 - -! input - integer, dimension(NGLLX,NGLLZ,nspec) :: ibool - -! output - integer, dimension(nspec) :: perm - -! global corner numbers that need to be created - integer, dimension(nglob) :: global_corner_number - - integer mn(nspec*NGNOD_QUADRANGLE),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 - -! 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 - -! perform inverse Cuthill-McKee (1969) permutation for mesh numbering - logical :: PERFORM_CUTHILL_MCKEE - -! -!----------------------------------------------------------------------- -! - 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_four_corners_only) - do i=1,nspec - istart = mp(i) - istop = mp(i+1) - 1 - enddo - - allocate(np(nglob_four_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_four_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_four_corners_only,nspec,count_only,total_size_ne) - do i=1,nglob_four_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_four_corners_only,& - count_only,total_size_ne,total_size_adj) - deallocate(adj) - allocate(adj(total_size_adj)) - count_only = .false. - call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_four_corners_only,& - count_only,total_size_ne,total_size_adj) - do i=1,nspec - istart = xadj(i) - istop = xadj(i+1) - 1 - number_of_neighbors = istop-istart+1 - if(number_of_neighbors < 1) stop 'error: your mesh seems to have at least one element not connected to any other' - if(number_of_neighbors > MAX_NUMBER_OF_NEIGHBORS) stop 'error: your mesh seems to have an unlikely high valence' - enddo - deallocate(ne,np) - -! call the Cuthill-McKee sorting algorithm - call cuthill_mckee(adj,xadj,perm,invperm,nspec,total_size_adj,limit) - deallocate(adj) - else -! create identity permutation in order to do nothing - do i=1,nspec - perm(i) = i - enddo - endif - - end subroutine get_perm_cuthill_mckee - -!======================================================================= -! -! 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_four_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_QUADRANGLE 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,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_QUADRANGLE),mp(nspec+1) - integer, intent(out) :: nglob_four_corners_only - - integer ninter,nsum,ispec,node,k,inumcorner,ix,iy - - 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_four_corners_only = 0 - global_corner_number(:) = -1 - - do ispec=1,nspec - - inumcorner = 0 - do iy = 1,NGLLZ,NGLLZ-1 - do ix = 1,NGLLX,NGLLX-1 - - inumcorner = inumcorner + 1 - if(inumcorner > NGNOD_QUADRANGLE) 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,ispec)) == -1) then - nglob_four_corners_only = nglob_four_corners_only + 1 - global_corner_number(ibool(ix,iy,ispec)) = nglob_four_corners_only - endif - - node = global_corner_number(ibool(ix,iy,ispec)) - do k=nsum,ninter-1 - if(node == mn(k)) goto 200 - enddo - - mn(ninter) = node - ninter = ninter + 1 - 200 continue - - 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_four_corners_only, & - nspec,count_only,total_size_ne) - -!----------------------------------------------------------------------- -! -! Forms the NE and NP arrays -! -! Input : -! ------- -! MN, MP, nspec -! nglob_four_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_four_corners_only,nspec - - integer, intent(in) :: mn(nspec*NGNOD_QUADRANGLE),mp(nspec+1) - - integer, intent(out) :: ne(total_size_ne),np(nglob_four_corners_only+1) - - integer nsum,inode,ispec,j - - nsum = 1 - np(1) = 1 - - do inode=1,nglob_four_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_four_corners_only,count_only,total_size_ne,total_size_adj) - -!----------------------------------------------------------------------- -! -! 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 - integer total_size_ne,total_size_adj - - integer nglob_four_corners_only - - integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel - - integer, intent(in) :: mn(nspec*NGNOD_QUADRANGLE),mp(nspec+1),ne(total_size_ne),np(nglob_four_corners_only+1) - - integer, intent(out) :: adj(total_size_adj),xadj(nspec+1) - - logical maskel(nspec) - integer countel(nspec) - - 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 -!! DK DK this below implemented by David Michea in 3D, but not true anymore in 2D: should be -!! DK DK two corners instead of three. But does not matter because FACE is always .false. -!! DK DK and therefore this part of the routine is currently never used. -!! DK DK Let me add a stop statement just in case. - stop 'FACE = .true. not implemented, check the above comment in the source code' -!! DK DK End of the stop statement added. - ! 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) - - 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 - 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 elements 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 permutation (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 permutation (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 - -! added a test to detect disconnected subsets in the mesh -! (in which case Cuthill-McKee fails and should be turned off) - if(root > nbnodes+1) stop 'error: root > nbnodes+1 in Cuthill-McKee' - if(root < 1) then - print *,'error: root < 1 in Cuthill-McKee; you probably have a mesh composed of' - print *,'two disconnected subsets of elements, in which case Cuthill-McKee fails and should be turned off.' - print *,'please set PERFORM_CUTHILL_MCKEE = .false. in constants.h and recompile.' - print *,'please also doublecheck that you indeed want to run two separate meshes simultaneously,' - print *,'which is extremely unusual (but formally not incorrect).' - stop 'fatal error in Cuthill-McKee' - endif - -! 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 - diff --git a/src/specfem2D/gmat01.f90 b/src/specfem2D/gmat01.f90 index 513cfabe3..3a87bf505 100644 --- a/src/specfem2D/gmat01.f90 +++ b/src/specfem2D/gmat01.f90 @@ -44,7 +44,7 @@ subroutine gmat01(density_array,porosity_array,tortuosity_array, & aniso_array,permeability,poroelastcoef, & - numat,myrank,ipass,QKappa_array,Qmu_array, & + numat,myrank,QKappa_array,Qmu_array, & freq0,Q0,f0,ATTENUATION_PORO_FLUID_PART) ! reads properties of a 2D isotropic or anisotropic linear elastic element @@ -52,7 +52,7 @@ subroutine gmat01(density_array,porosity_array,tortuosity_array, & implicit none include "constants.h" - integer :: numat,myrank,ipass + integer :: numat,myrank double precision :: density_array(2,numat),poroelastcoef(4,3,numat),porosity_array(numat) double precision :: aniso_array(9,numat),tortuosity_array(numat),permeability(3,numat) double precision :: QKappa_array(numat),Qmu_array(numat) @@ -86,7 +86,7 @@ subroutine gmat01(density_array,porosity_array,tortuosity_array, & QKappa_array(:) = 9999. Qmu_array(:) = 9999. - if(myrank == 0 .and. ipass == 1) write(IOUT,100) numat + if(myrank == 0) write(IOUT,100) numat read(IIN,"(a80)") datlin read(IIN,"(a80)") datlin @@ -277,7 +277,7 @@ subroutine gmat01(density_array,porosity_array,tortuosity_array, & ! !---- check what has been read ! - if(myrank == 0 .and. ipass == 1) then + if(myrank == 0) then if(indic == 1) then ! material can be acoustic (fluid) or elastic (solid) if(poroelastcoef(2,1,n) > TINYVAL) then ! elastic diff --git a/src/specfem2D/initialize_simulation.F90 b/src/specfem2D/initialize_simulation.F90 index d3f3ab0de..9038775d1 100644 --- a/src/specfem2D/initialize_simulation.F90 +++ b/src/specfem2D/initialize_simulation.F90 @@ -44,8 +44,7 @@ !======================================================================== - subroutine initialize_simulation(nproc,myrank,NUMBER_OF_PASSES, & - ninterface_acoustic,ninterface_elastic,ninterface_poroelastic) + subroutine initialize_simulation(nproc,myrank,ninterface_acoustic,ninterface_elastic,ninterface_poroelastic) #ifdef USE_MPI use :: mpi @@ -53,7 +52,7 @@ subroutine initialize_simulation(nproc,myrank,NUMBER_OF_PASSES, & implicit none include "constants.h" - integer :: nproc,myrank,NUMBER_OF_PASSES + integer :: nproc,myrank integer :: ninterface_acoustic, ninterface_elastic,ninterface_poroelastic ! local parameters @@ -71,18 +70,9 @@ subroutine initialize_simulation(nproc,myrank,NUMBER_OF_PASSES, & call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ier) call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier) if( ier /= 0 ) call exit_MPI('error MPI initialization') - - ! this is only used in the case of MPI because it distinguishes between inner and outer element - ! in the MPI partitions, which is meaningless in the serial case - if(FURTHER_REDUCE_CACHE_MISSES) then - NUMBER_OF_PASSES = 2 - else - NUMBER_OF_PASSES = 1 - endif #else nproc = 1 myrank = 0 - NUMBER_OF_PASSES = 1 ! assign dummy value for now, will be changed later if needed #endif ninterface_acoustic = 0 diff --git a/src/specfem2D/locate_receivers.F90 b/src/specfem2D/locate_receivers.F90 index 5c7068b66..31e046eca 100644 --- a/src/specfem2D/locate_receivers.F90 +++ b/src/specfem2D/locate_receivers.F90 @@ -51,7 +51,7 @@ subroutine locate_receivers(ibool,coord,nspec,nglob,xigll,zigll, & st_xval,st_zval,ispec_selected_rec, & xi_receiver,gamma_receiver,station_name,network_name, & x_source,z_source, & - coorg,knods,ngnod,npgeo,ipass, & + coorg,knods,ngnod,npgeo, & x_final_receiver, z_final_receiver) #ifdef USE_MPI @@ -62,7 +62,7 @@ subroutine locate_receivers(ibool,coord,nspec,nglob,xigll,zigll, & include "constants.h" - integer nrec,nspec,nglob,ngnod,npgeo,ipass + integer nrec,nspec,nglob,ngnod,npgeo integer, intent(in) :: nproc, myrank integer knods(ngnod,nspec) @@ -117,7 +117,7 @@ subroutine locate_receivers(ibool,coord,nspec,nglob,xigll,zigll, & ! ************** - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then write(IOUT,*) write(IOUT,*) '********************' write(IOUT,*) ' locating receivers' @@ -271,7 +271,7 @@ subroutine locate_receivers(ibool,coord,nspec,nglob,xigll,zigll, & endif enddo - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then do irec = 1, nrec write(IOUT,*) diff --git a/src/specfem2D/locate_source_force.F90 b/src/specfem2D/locate_source_force.F90 index 3c1a27cef..bca2e6c28 100644 --- a/src/specfem2D/locate_source_force.F90 +++ b/src/specfem2D/locate_source_force.F90 @@ -48,7 +48,7 @@ subroutine locate_source_force(ibool,coord,nspec,nglob,xigll,zigll,x_source,z_source, & ispec_selected_source,is_proc_source,nb_proc_source,nproc,myrank, & - xi_source,gamma_source,coorg,knods,ngnod,npgeo,ipass,iglob_source) + xi_source,gamma_source,coorg,knods,ngnod,npgeo,iglob_source) #ifdef USE_MPI use :: mpi @@ -58,7 +58,7 @@ subroutine locate_source_force(ibool,coord,nspec,nglob,xigll,zigll,x_source,z_so include "constants.h" - integer nspec,nglob,ngnod,npgeo,ipass + integer nspec,nglob,ngnod,npgeo integer knods(ngnod,nspec) double precision coorg(NDIM,npgeo) @@ -94,7 +94,7 @@ subroutine locate_source_force(ibool,coord,nspec,nglob,xigll,zigll,x_source,z_so ! ************** - if ((myrank == 0 .or. nproc == 1) .and. ipass == 1) then + if (myrank == 0 .or. nproc == 1) then write(IOUT,*) write(IOUT,*) '*******************************' write(IOUT,*) ' locating force source' @@ -229,7 +229,7 @@ subroutine locate_source_force(ibool,coord,nspec,nglob,xigll,zigll,x_source,z_so ! compute final distance between asked and found final_distance = sqrt((x_source-x)**2 + (z_source-z)**2) - if (is_proc_source == 1 .and. ipass == 1) then + if (is_proc_source == 1) then write(IOUT,*) write(IOUT,*) 'Force source:' diff --git a/src/specfem2D/locate_source_moment_tensor.F90 b/src/specfem2D/locate_source_moment_tensor.F90 index 3bf8aa61d..b392ed6ba 100644 --- a/src/specfem2D/locate_source_moment_tensor.F90 +++ b/src/specfem2D/locate_source_moment_tensor.F90 @@ -49,7 +49,7 @@ subroutine locate_source_moment_tensor(ibool,coord,nspec,nglob, & xigll,zigll,x_source,z_source, & ispec_selected_source,is_proc_source,nb_proc_source,nproc,myrank, & - xi_source,gamma_source,coorg,knods,ngnod,npgeo,ipass) + xi_source,gamma_source,coorg,knods,ngnod,npgeo) #ifdef USE_MPI use :: mpi @@ -59,7 +59,7 @@ subroutine locate_source_moment_tensor(ibool,coord,nspec,nglob, & include "constants.h" - integer nspec,nglob,ngnod,npgeo,ipass + integer nspec,nglob,ngnod,npgeo integer knods(ngnod,nspec) double precision coorg(NDIM,npgeo) @@ -92,10 +92,9 @@ subroutine locate_source_moment_tensor(ibool,coord,nspec,nglob, & integer :: ierror #endif - - ! ************** - if ((myrank == 0 .or. nproc == 1) .and. ipass == 1) then + + if (myrank == 0 .or. nproc == 1) then write(IOUT,*) write(IOUT,*) '*******************************' write(IOUT,*) ' locating moment-tensor source' @@ -228,7 +227,7 @@ subroutine locate_source_moment_tensor(ibool,coord,nspec,nglob, & ! compute final distance between asked and found final_distance = sqrt((x_source-x)**2 + (z_source-z)**2) - if (is_proc_source == 1 .and. ipass == 1) then + if (is_proc_source == 1) then write(IOUT,*) write(IOUT,*) 'Moment-tensor source:' diff --git a/src/specfem2D/read_databases.F90 b/src/specfem2D/read_databases.F90 index 4b28f3714..af2335a00 100644 --- a/src/specfem2D/read_databases.F90 +++ b/src/specfem2D/read_databases.F90 @@ -43,7 +43,7 @@ ! !======================================================================== - subroutine read_databases_init(myrank,ipass, & + subroutine read_databases_init(myrank, & simulation_title,SIMULATION_TYPE,NOISE_TOMOGRAPHY,SAVE_FORWARD,npgeo,nproc, & output_grid_Gnuplot,interpol,NSTEP_BETWEEN_OUTPUT_INFO,NSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP_BETWEEN_OUTPUT_IMAGES, & PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,NELEM_PML_THICKNESS, & @@ -57,7 +57,7 @@ subroutine read_databases_init(myrank,ipass, & save_binary_seismograms_single,save_binary_seismograms_double, & DRAW_SOURCES_AND_RECEIVERS,Q0,freq0,p_sv,NSTEP,deltat,NSOURCES, & factor_subsample_image,USE_SNAPSHOT_NUMBER_IN_FILENAME,DRAW_WATER_IN_BLUE,US_LETTER, & - POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0,time_stepping_scheme,& + POWER_DISPLAY_COLOR,SU_FORMAT,USER_T0,time_stepping_scheme,& ADD_SPRING_TO_STACEY,ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL,& read_external_mesh,save_ASCII_kernels) @@ -66,7 +66,7 @@ subroutine read_databases_init(myrank,ipass, & implicit none include "constants.h" - integer :: myrank,ipass + integer :: myrank character(len=60) simulation_title integer :: SIMULATION_TYPE,NOISE_TOMOGRAPHY,npgeo,nproc integer :: colors,numbers,subsamp_postscript,seismotype,imagetype_postscript @@ -104,9 +104,6 @@ subroutine read_databases_init(myrank,ipass, & ! non linear display to enhance small amplitudes in color images double precision :: POWER_DISPLAY_COLOR -! perform inverse Cuthill-McKee (1969) permutation for mesh numbering - logical :: PERFORM_CUTHILL_MCKEE - ! output seismograms in Seismic Unix format (adjoint traces will be read in the same format) logical :: SU_FORMAT @@ -154,9 +151,9 @@ subroutine read_databases_init(myrank,ipass, & read(IIN,"(a50)") simulation_title !---- print the date, time and start-up banner - if (myrank == 0 .and. ipass == 1) call datim(simulation_title) + if (myrank == 0) call datim(simulation_title) - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then write(IOUT,*) write(IOUT,*) write(IOUT,*) '*********************' @@ -279,9 +276,6 @@ subroutine read_databases_init(myrank,ipass, & read(IIN,"(a80)") datlin read(IIN,*) POWER_DISPLAY_COLOR - read(IIN,"(a80)") datlin - read(IIN,*) PERFORM_CUTHILL_MCKEE - read(IIN,"(a80)") datlin read(IIN,*) SU_FORMAT @@ -304,7 +298,7 @@ subroutine read_databases_init(myrank,ipass, & read(IIN,*) PERIODIC_DETECT_TOL !---- check parameters read - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then write(IOUT,200) npgeo,NDIM write(IOUT,600) NSTEP_BETWEEN_OUTPUT_INFO,colors,numbers write(IOUT,700) seismotype,anglerec @@ -317,7 +311,7 @@ subroutine read_databases_init(myrank,ipass, & !---- read time step read(IIN,"(a80)") datlin read(IIN,*) NSTEP,deltat - if (myrank == 0 .and. ipass == 1) write(IOUT,703) NSTEP,deltat,NSTEP*deltat + if (myrank == 0) write(IOUT,703) NSTEP,deltat,NSTEP*deltat if(SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. & (ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART)) then @@ -459,7 +453,7 @@ end subroutine read_databases_atten !------------------------------------------------------------------------------------------------- ! - subroutine read_databases_coorg_elem(myrank,ipass,npgeo,coorg,numat,ngnod,nspec, & + subroutine read_databases_coorg_elem(myrank,npgeo,coorg,numat,ngnod,nspec, & pointsdisp,plot_lowerleft_corner_only, & nelemabs,nelem_acoustic_surface, & num_fluid_solid_edges,num_fluid_poro_edges, & @@ -470,7 +464,7 @@ subroutine read_databases_coorg_elem(myrank,ipass,npgeo,coorg,numat,ngnod,nspec, implicit none include "constants.h" - integer :: myrank,ipass,npgeo + integer :: myrank,npgeo double precision, dimension(NDIM,npgeo) :: coorg integer :: numat,ngnod,nspec @@ -516,7 +510,7 @@ subroutine read_databases_coorg_elem(myrank,ipass,npgeo,coorg,numat,ngnod,nspec, num_solid_poro_edges,nnodes_tangential_curve !---- print element group main parameters - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0 ) then write(IOUT,107) write(IOUT,207) nspec,ngnod,NGLLX,NGLLZ,NGLLX*NGLLZ,pointsdisp,numat endif @@ -539,21 +533,18 @@ end subroutine read_databases_coorg_elem !------------------------------------------------------------------------------------------------- ! - subroutine read_databases_mato(ipass,nspec,ngnod,kmato,knods, & - perm,antecedent_list,region_CPML) + subroutine read_databases_mato(nspec,ngnod,kmato,knods,region_CPML) ! reads spectral macrobloc data implicit none include "constants.h" - integer :: ipass,ngnod,nspec + integer :: ngnod,nspec integer, dimension(nspec) :: kmato integer, dimension(nspec) :: region_CPML integer, dimension(ngnod,nspec) :: knods - integer, dimension(nspec) :: perm,antecedent_list - ! local parameters integer :: n,k,ispec,kmato_read,pml_read integer, dimension(:), allocatable :: knods_read @@ -573,19 +564,11 @@ subroutine read_databases_mato(ipass,nspec,ngnod,kmato,knods, & do ispec = 1,nspec ! format: #element_id #material_id #node_id1 #node_id2 #... read(IIN,*) n,kmato_read,(knods_read(k), k=1,ngnod),pml_read - if(ipass == 1) then ! material association kmato(n) = kmato_read region_CPML(n) = pml_read ! element control node indices knods(:,n)= knods_read(:) - else if(ipass == 2) then - kmato(perm(antecedent_list(n))) = kmato_read - region_CPML(perm(antecedent_list(n))) = pml_read - knods(:,perm(antecedent_list(n)))= knods_read(:) - else - call exit_MPI('error: maximum is 2 passes') - endif enddo deallocate(knods_read) @@ -617,22 +600,17 @@ end subroutine read_databases_ninterface !------------------------------------------------------------------------------------------------- ! - subroutine read_databases_interfaces(ipass,ninterface,nspec,max_interface_size, & - my_neighbours,my_nelmnts_neighbours,my_interfaces, & - perm,antecedent_list) + subroutine read_databases_interfaces(ninterface,max_interface_size,my_neighbours,my_nelmnts_neighbours,my_interfaces) ! reads in interfaces implicit none include "constants.h" - integer :: ipass,nspec integer :: ninterface,max_interface_size integer, dimension(ninterface) :: my_neighbours,my_nelmnts_neighbours integer, dimension(4,max_interface_size,ninterface) :: my_interfaces - integer, dimension(nspec) :: perm,antecedent_list - ! local parameters integer :: num_interface,ie,my_interfaces_read @@ -659,13 +637,7 @@ subroutine read_databases_interfaces(ipass,ninterface,nspec,max_interface_size, read(IIN,*) my_interfaces_read, my_interfaces(2,ie,num_interface), & my_interfaces(3,ie,num_interface), my_interfaces(4,ie,num_interface) - if(ipass == 1) then - my_interfaces(1,ie,num_interface) = my_interfaces_read - else if(ipass == 2) then - my_interfaces(1,ie,num_interface) = perm(antecedent_list(my_interfaces_read)) - else - call exit_MPI('error: maximum number of passes is 2') - endif + my_interfaces(1,ie,num_interface) = my_interfaces_read enddo enddo @@ -677,10 +649,10 @@ end subroutine read_databases_interfaces !------------------------------------------------------------------------------------------------- ! - subroutine read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, & + subroutine read_databases_absorbing(myrank,nelemabs,nspec,anyabs, & ibegin_edge1,iend_edge1,ibegin_edge2,iend_edge2, & ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4, & - numabs,codeabs,typeabs,perm,antecedent_list, & + numabs,codeabs,typeabs, & nspec_left,nspec_right,nspec_bottom,nspec_top, & ib_right,ib_left,ib_bottom,ib_top,PML_BOUNDARY_CONDITIONS) @@ -693,14 +665,13 @@ subroutine read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, & implicit none include "constants.h" - integer :: myrank,ipass,nspec + integer :: myrank,nspec integer :: nelemabs integer, dimension(nelemabs) :: numabs,ibegin_edge1,iend_edge1, & ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4,ibegin_edge2,iend_edge2 logical, dimension(4,nelemabs) :: codeabs integer, dimension(nelemabs) :: typeabs logical :: anyabs,PML_BOUNDARY_CONDITIONS - integer, dimension(nspec) :: perm,antecedent_list integer :: nspec_left,nspec_right,nspec_bottom,nspec_top integer, dimension(nelemabs) :: ib_right,ib_left,ib_bottom,ib_top @@ -761,13 +732,7 @@ subroutine read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, & if(numabsread < 1 .or. numabsread > nspec) & call exit_MPI('Wrong absorbing element number') - if(ipass == 1) then - numabs(inum) = numabsread - else if(ipass == 2) then - numabs(inum) = perm(antecedent_list(numabsread)) - else - call exit_MPI('error: maximum number of passes is 2') - endif + numabs(inum) = numabsread codeabs(IEDGE1,inum) = codeabsread(1) codeabs(IEDGE2,inum) = codeabsread(2) @@ -811,35 +776,29 @@ subroutine read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, & else ! if this MPI slice has no absorbing element at all - if (ipass == 1) then nelemabs = 0 nspec_left = 0 nspec_right = 0 nspec_bottom = 0 nspec_top = 0 - endif endif #ifdef USE_MPI - if (ipass == 1) then call MPI_REDUCE(nelemabs, nelemabs_tot, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier) call MPI_REDUCE(nspec_left, nspec_left_tot, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier) call MPI_REDUCE(nspec_right, nspec_right_tot, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier) call MPI_REDUCE(nspec_bottom, nspec_bottom_tot, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier) call MPI_REDUCE(nspec_top, nspec_top_tot, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ier) - endif #else - if (ipass == 1) then nelemabs_tot = nelemabs nspec_left_tot = nspec_left nspec_right_tot = nspec_right nspec_bottom_tot = nspec_bottom nspec_top_tot = nspec_top - endif #endif - if (myrank == 0 .and. ipass == 1 .and. .not. PML_BOUNDARY_CONDITIONS) then + if (myrank == 0 .and. .not. PML_BOUNDARY_CONDITIONS) then write(IOUT,*) write(IOUT,*) 'Number of absorbing elements: ',nelemabs_tot write(IOUT,*) ' nspec_left = ',nspec_left_tot @@ -855,21 +814,17 @@ end subroutine read_databases_absorbing !------------------------------------------------------------------------------------------------- ! - subroutine read_databases_free_surf(ipass,nelem_acoustic_surface,nspec, & - acoustic_edges,perm,antecedent_list,any_acoustic_edges) + subroutine read_databases_free_surf(nelem_acoustic_surface,acoustic_edges,any_acoustic_edges) ! reads acoustic free surface data implicit none include "constants.h" - integer :: ipass,nspec integer :: nelem_acoustic_surface integer, dimension(4,nelem_acoustic_surface) :: acoustic_edges logical :: any_acoustic_edges - integer, dimension(nspec) :: perm,antecedent_list - ! local parameters integer :: inum,acoustic_edges_read character(len=80) :: datlin @@ -885,13 +840,7 @@ subroutine read_databases_free_surf(ipass,nelem_acoustic_surface,nspec, & read(IIN,*) acoustic_edges_read, acoustic_edges(2,inum), acoustic_edges(3,inum), & acoustic_edges(4,inum) - if(ipass == 1) then - acoustic_edges(1,inum) = acoustic_edges_read - else if(ipass == 2) then - acoustic_edges(1,inum) = perm(antecedent_list(acoustic_edges_read)) - else - call exit_MPI('error: maximum number of passes is 2') - endif + acoustic_edges(1,inum) = acoustic_edges_read enddo @@ -903,13 +852,12 @@ end subroutine read_databases_free_surf !------------------------------------------------------------------------------------------------- ! - subroutine read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_solid_edges, & + subroutine read_databases_coupled(num_fluid_solid_edges,any_fluid_solid_edges, & fluid_solid_acoustic_ispec,fluid_solid_elastic_ispec, & num_fluid_poro_edges,any_fluid_poro_edges, & fluid_poro_acoustic_ispec,fluid_poro_poroelastic_ispec, & num_solid_poro_edges,any_solid_poro_edges, & - solid_poro_elastic_ispec,solid_poro_poroelastic_ispec, & - perm,antecedent_list) + solid_poro_elastic_ispec,solid_poro_poroelastic_ispec) ! reads acoustic elastic coupled edges ! reads acoustic poroelastic coupled edges @@ -918,8 +866,6 @@ subroutine read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_so implicit none include "constants.h" - integer :: ipass,nspec - integer :: num_fluid_solid_edges logical :: any_fluid_solid_edges integer, dimension(num_fluid_solid_edges) :: fluid_solid_acoustic_ispec,fluid_solid_elastic_ispec @@ -932,8 +878,6 @@ subroutine read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_so logical :: any_solid_poro_edges integer, dimension(num_solid_poro_edges) :: solid_poro_elastic_ispec,solid_poro_poroelastic_ispec - integer, dimension(nspec) :: perm,antecedent_list - ! local parameters integer :: inum integer :: fluid_solid_acoustic_ispec_read,fluid_solid_elastic_ispec_read, & @@ -956,15 +900,8 @@ subroutine read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_so do inum = 1, num_fluid_solid_edges read(IIN,*) fluid_solid_acoustic_ispec_read,fluid_solid_elastic_ispec_read - if(ipass == 1) then - fluid_solid_acoustic_ispec(inum) = fluid_solid_acoustic_ispec_read - fluid_solid_elastic_ispec(inum) = fluid_solid_elastic_ispec_read - else if(ipass == 2) then - fluid_solid_acoustic_ispec(inum) = perm(antecedent_list(fluid_solid_acoustic_ispec_read)) - fluid_solid_elastic_ispec(inum) = perm(antecedent_list(fluid_solid_elastic_ispec_read)) - else - call exit_MPI('error: maximum number of passes is 2') - endif + fluid_solid_acoustic_ispec(inum) = fluid_solid_acoustic_ispec_read + fluid_solid_elastic_ispec(inum) = fluid_solid_elastic_ispec_read enddo endif @@ -975,15 +912,8 @@ subroutine read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_so do inum = 1, num_fluid_poro_edges read(IIN,*) fluid_poro_acoustic_ispec_read,fluid_poro_poro_ispec_read - if(ipass == 1) then - fluid_poro_acoustic_ispec(inum) = fluid_poro_acoustic_ispec_read - fluid_poro_poroelastic_ispec(inum) = fluid_poro_poro_ispec_read - else if(ipass == 2) then - fluid_poro_acoustic_ispec(inum) = perm(antecedent_list(fluid_poro_acoustic_ispec_read)) - fluid_poro_poroelastic_ispec(inum) = perm(antecedent_list(fluid_poro_poro_ispec_read)) - else - call exit_MPI('error: maximum number of passes is 2') - endif + fluid_poro_acoustic_ispec(inum) = fluid_poro_acoustic_ispec_read + fluid_poro_poroelastic_ispec(inum) = fluid_poro_poro_ispec_read enddo endif @@ -994,15 +924,8 @@ subroutine read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_so do inum = 1, num_solid_poro_edges read(IIN,*) solid_poro_poro_ispec_read,solid_poro_elastic_ispec_read - if(ipass == 1) then - solid_poro_elastic_ispec(inum) = solid_poro_elastic_ispec_read - solid_poro_poroelastic_ispec(inum) = solid_poro_poro_ispec_read - else if(ipass == 2) then - solid_poro_elastic_ispec(inum) = perm(antecedent_list(solid_poro_elastic_ispec_read)) - solid_poro_poroelastic_ispec(inum) = perm(antecedent_list(solid_poro_poro_ispec_read)) - else - call exit_MPI('error: maximum number of passes is 2') - endif + solid_poro_elastic_ispec(inum) = solid_poro_elastic_ispec_read + solid_poro_poroelastic_ispec(inum) = solid_poro_poro_ispec_read enddo endif @@ -1014,7 +937,7 @@ end subroutine read_databases_coupled subroutine read_databases_final(nnodes_tangential_curve,nodes_tangential_curve, & force_normal_to_surface,rec_normal_to_surface, & - any_tangential_curve ) + any_tangential_curve) ! reads tangential detection curve ! and closes Database file diff --git a/src/specfem2D/set_sources.f90 b/src/specfem2D/set_sources.f90 index b1c365717..41bb0084d 100644 --- a/src/specfem2D/set_sources.f90 +++ b/src/specfem2D/set_sources.f90 @@ -46,7 +46,7 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & x_source,z_source,Mxx,Mzz,Mxz,f0,tshift_src,factor,anglesource,aval, & - t0,initialfield,ipass,deltat,USER_T0) + t0,initialfield,deltat,USER_T0) ! gets source parameters @@ -61,7 +61,6 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & double precision, dimension(NSOURCES) :: aval double precision :: t0 double precision :: deltat - integer :: ipass logical :: initialfield ! use this t0 as earliest starting time rather than the automatically calculated one @@ -80,13 +79,13 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & ! checks source type if(.not. initialfield) then if (source_type(i_source) == 1) then - if ( myrank == 0 .and. ipass == 1 ) then + if ( myrank == 0 ) then ! user output write(IOUT,212) x_source(i_source),z_source(i_source),f0(i_source),tshift_src(i_source), & factor(i_source),anglesource(i_source) endif else if(source_type(i_source) == 2) then - if ( myrank == 0 .and. ipass == 1 ) then + if ( myrank == 0 ) then ! user output write(IOUT,222) x_source(i_source),z_source(i_source),f0(i_source),tshift_src(i_source), & factor(i_source),Mxx(i_source),Mzz(i_source),Mxz(i_source) @@ -147,7 +146,7 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & ! time 0 on time axis will correspond to given origin time ! notifies user - if( myrank == 0 .and. ipass == 1) then + if( myrank == 0 ) then write(IOUT,*) write(IOUT,*) ' using USER_T0 . . . . . . . . . = ',USER_T0 write(IOUT,*) ' original t0 . . . . . . . . . = ',t0 @@ -164,7 +163,7 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & t0 = USER_T0 ! notifies user - if( myrank == 0 .and. ipass == 1) then + if( myrank == 0 ) then write(IOUT,*) ' fix new simulation start time . = ', - t0 endif @@ -177,19 +176,19 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & tshift_src(i_source) = t0_source(i_source) - 1.20d0 * hdur(i_source) endif ! user output - if( myrank == 0 .and. ipass == 1) then + if( myrank == 0 ) then write(IOUT,*) ' source ',i_source,'uses tshift = ',tshift_src(i_source) endif enddo ! user output - if( myrank == 0 .and. ipass == 1) then + if( myrank == 0 ) then write(IOUT,*) endif else ! start time needs to be at least t0 for numerical stability ! notifies user - if( myrank == 0 .and. ipass == 1) then + if( myrank == 0 ) then write(IOUT,*) 'error: USER_T0 is too small' write(IOUT,*) ' must make one of three adjustements:' write(IOUT,*) ' - increase USER_T0 to be at least: ',t0 @@ -199,7 +198,7 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & call exit_MPI('error USER_T0 is set but too small') endif else if( USER_T0 < 0.d0 ) then - if( myrank == 0 .and. ipass == 1 ) then + if( myrank == 0 ) then write(IOUT,*) 'error: USER_T0 is negative, must be set zero or positive!' endif call exit_MPI('error negative USER_T0 parameter in constants.h') @@ -215,7 +214,7 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & if(time_function_type(i_source) /= 4 .and. time_function_type(i_source) /= 5) then ! user output - if( myrank == 0 .and. ipass == 1 ) then + if( myrank == 0 ) then write(IOUT,*) ' Onset time. . . . . . = ',t0+tshift_src(i_source) write(IOUT,*) ' Fundamental period. . = ',1.d0/f0(i_source) write(IOUT,*) ' Fundamental frequency = ',f0(i_source) @@ -225,7 +224,7 @@ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, & if( t0+tshift_src(i_source) <= 1.d0/f0(i_source)) then call exit_MPI('Onset time too small') else - if( myrank == 0 .and. ipass == 1 ) then + if( myrank == 0 ) then write(IOUT,*) ' --> onset time ok' endif endif diff --git a/src/specfem2D/setup_sources_receivers.F90 b/src/specfem2D/setup_sources_receivers.F90 index 640bd32d6..e960eadb4 100644 --- a/src/specfem2D/setup_sources_receivers.F90 +++ b/src/specfem2D/setup_sources_receivers.F90 @@ -43,11 +43,11 @@ ! !======================================================================== - subroutine setup_sources_receivers(NSOURCES,initialfield,source_type,& + subroutine setup_sources_receivers(NSOURCES,initialfield,source_type, & coord,ibool,nglob,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, & x_source,z_source,ispec_selected_source,ispec_selected_rec, & - is_proc_source,nb_proc_source,ipass,& - sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,& + is_proc_source,nb_proc_source, & + sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo, & nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod, & nrec,nrecloc,recloc,which_proc_receiver,st_xval,st_zval, & xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver,iglob_source) @@ -58,7 +58,7 @@ subroutine setup_sources_receivers(NSOURCES,initialfield,source_type,& logical :: initialfield integer :: NSOURCES - integer :: npgeo,ngnod,myrank,ipass,nproc + integer :: npgeo,ngnod,myrank,nproc integer :: nglob,nspec,nelem_acoustic_surface ! Gauss-Lobatto-Legendre points @@ -104,8 +104,8 @@ subroutine setup_sources_receivers(NSOURCES,initialfield,source_type,& ! collocated force source call locate_source_force(ibool,coord,nspec,nglob,xigll,zigll,x_source(i_source),z_source(i_source), & - ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),& - nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass,& + ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source), & + nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo, & iglob_source(i_source)) ! check that acoustic source is not exactly on the free surface because pressure is zero there @@ -144,8 +144,8 @@ subroutine setup_sources_receivers(NSOURCES,initialfield,source_type,& else if(source_type(i_source) == 2) then ! moment-tensor source call locate_source_moment_tensor(ibool,coord,nspec,nglob,xigll,zigll,x_source(i_source),z_source(i_source), & - ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),& - nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass) + ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source), & + nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo) ! compute source array for moment-tensor source call compute_arrays_source(ispec_selected_source(i_source),xi_source(i_source),gamma_source(i_source),& @@ -166,14 +166,14 @@ subroutine setup_sources_receivers(NSOURCES,initialfield,source_type,& st_xval,st_zval,ispec_selected_rec, & xi_receiver,gamma_receiver,station_name,network_name, & x_source(1),z_source(1), & - coorg,knods,ngnod,npgeo,ipass, & + coorg,knods,ngnod,npgeo, & x_final_receiver,z_final_receiver) !! DK DK this below not supported in the case of MPI yet, we should do a MPI_GATHER() of the values !! DK DK and use "if(myrank == which_proc_receiver(irec)) then" to display the right sources !! DK DK and receivers carried by each mesh slice, and not fictitious values coming from other slices #ifndef USE_MPI - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then ! write actual source locations to file ! note that these may differ from input values, especially if source_surf = .true. in SOURCE diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90 index cc1aeaf06..09df9654f 100644 --- a/src/specfem2D/specfem2D.F90 +++ b/src/specfem2D/specfem2D.F90 @@ -403,9 +403,6 @@ program specfem2D ! non linear display to enhance small amplitudes in color images double precision :: POWER_DISPLAY_COLOR -! perform inverse Cuthill-McKee (1969) permutation for mesh numbering - logical :: PERFORM_CUTHILL_MCKEE - ! output seismograms in Seismic Unix format (adjoint traces will be read in the same format) logical :: SU_FORMAT @@ -823,10 +820,9 @@ program specfem2D integer count_left,count_right,count_bottom logical :: over_critical_angle -! further reduce cache misses inner/outer in two passes in the case of an MPI simulation - integer :: ipass,ispec_inner,ispec_outer,NUMBER_OF_PASSES +! inner/outer elements in the case of an MPI simulation + integer :: ispec_inner,ispec_outer integer :: nglob_outer,nglob_inner - integer, dimension(:), allocatable :: perm,antecedent_list,check_perm ! arrays for plotpost integer :: d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, & @@ -1061,15 +1057,11 @@ program specfem2D ! force Flush-To-Zero if available to avoid very slow Gradual Underflow trapping call force_ftz() - call initialize_simulation(nproc,myrank,NUMBER_OF_PASSES, & - ninterface_acoustic,ninterface_elastic,ninterface_poroelastic) + call initialize_simulation(nproc,myrank,ninterface_acoustic,ninterface_elastic,ninterface_poroelastic) if(nproc < 1) stop 'should have nproc >= 1' ! starts reading in Database file - ! it is necessary to duplicate this call before the loop on ipass = 1,NUMBER_OF_PASSES - ! because we need to read the value of logical flag PERFORM_CUTHILL_MCKEE - ipass = 1 - call read_databases_init(myrank,ipass, & + call read_databases_init(myrank, & simulation_title,SIMULATION_TYPE,NOISE_TOMOGRAPHY,SAVE_FORWARD,npgeo,nproc_read_from_database, & output_grid_Gnuplot,interpol,NSTEP_BETWEEN_OUTPUT_INFO,NSTEP_BETWEEN_OUTPUT_SEISMOS, & NSTEP_BETWEEN_OUTPUT_IMAGES,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,NELEM_PML_THICKNESS, & @@ -1084,7 +1076,7 @@ program specfem2D save_binary_seismograms_single,save_binary_seismograms_double,DRAW_SOURCES_AND_RECEIVERS, & Q0,freq0,p_sv,NSTEP,deltat,NSOURCES, & factor_subsample_image,USE_SNAPSHOT_NUMBER_IN_FILENAME,DRAW_WATER_IN_BLUE,US_LETTER, & - POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0, time_stepping_scheme, & + POWER_DISPLAY_COLOR,SU_FORMAT,USER_T0, time_stepping_scheme, & ADD_SPRING_TO_STACEY,ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL,& read_external_mesh,save_ASCII_kernels) @@ -1097,42 +1089,9 @@ program specfem2D npgeo_ori = npgeo if(ADD_A_SMALL_CRACK_IN_THE_MEDIUM) npgeo = npgeo + NB_POINTS_TO_ADD_TO_NPGEO -#ifdef USE_MPI - if(PERFORM_CUTHILL_MCKEE) then - NUMBER_OF_PASSES = 2 - else - NUMBER_OF_PASSES = 1 - endif -#endif - - ! reduction of cache misses inner/outer in two passes - do ipass = 1,NUMBER_OF_PASSES - - ! starts reading in Database file - if(ipass > 1) & - call read_databases_init(myrank,ipass, & - simulation_title,SIMULATION_TYPE,NOISE_TOMOGRAPHY,SAVE_FORWARD,npgeo,nproc_read_from_database, & - output_grid_Gnuplot,interpol,NSTEP_BETWEEN_OUTPUT_INFO,NSTEP_BETWEEN_OUTPUT_SEISMOS, & - NSTEP_BETWEEN_OUTPUT_IMAGES,PML_BOUNDARY_CONDITIONS,NELEM_PML_THICKNESS, & - NSTEP_BETWEEN_OUTPUT_WAVE_DUMPS,subsamp_seismos, & - imagetype_JPEG,imagetype_wavefield_dumps, & - output_postscript_snapshot,output_color_image,colors,numbers, & - meshvect,modelvect,boundvect,cutsnaps,subsamp_postscript,sizemax_arrows, & - anglerec,initialfield,add_Bielak_conditions, & - seismotype,imagetype_postscript,assign_external_model,READ_EXTERNAL_SEP_FILE, & - output_grid_ASCII,output_energy,output_wavefield_dumps,use_binary_for_wavefield_dumps, & - ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,save_ASCII_seismograms, & - save_binary_seismograms_single,save_binary_seismograms_double,DRAW_SOURCES_AND_RECEIVERS, & - Q0,freq0,p_sv,NSTEP,deltat,NSOURCES, & - factor_subsample_image,USE_SNAPSHOT_NUMBER_IN_FILENAME,DRAW_WATER_IN_BLUE,US_LETTER, & - POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0, time_stepping_scheme, & - ADD_SPRING_TO_STACEY,ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL,& - read_external_mesh) - ! !--- source information ! - if(ipass == 1) then allocate( source_type(NSOURCES) ) allocate( time_function_type(NSOURCES) ) allocate( x_source(NSOURCES) ) @@ -1155,7 +1114,6 @@ program specfem2D allocate( is_proc_source(NSOURCES) ) allocate( nb_proc_source(NSOURCES) ) allocate( sourcearray(NSOURCES,NDIM,NGLLX,NGLLZ) ) - endif ! reads in source infos call read_databases_sources(NSOURCES,source_type,time_function_type, & @@ -1164,7 +1122,7 @@ program specfem2D ! sets source parameters call set_sources(myrank,NSOURCES,source_type,time_function_type, & x_source,z_source,Mxx,Mzz,Mxz,f0,tshift_src,factor,anglesource,aval, & - t0,initialfield,ipass,deltat,USER_T0) + t0,initialfield,deltat,USER_T0) !---- define time stepping scheme if(time_stepping_scheme == 1)then @@ -1184,20 +1142,19 @@ program specfem2D endif !---- read the spectral macrobloc nodal coordinates - if(ipass == 1) allocate(coorg(NDIM,npgeo)) + allocate(coorg(NDIM,npgeo)) ! reads the spectral macrobloc nodal coordinates ! and basic properties of the spectral elements - !! DK DK call read_databases_coorg_elem(myrank,ipass,npgeo,coorg,numat,ngnod,nspec, & + !! DK DK call read_databases_coorg_elem(myrank,npgeo,coorg,numat,ngnod,nspec, & !! DK DK Dec 2011: added a crack manually - call read_databases_coorg_elem(myrank,ipass,npgeo_ori,coorg,numat,ngnod,nspec, & + call read_databases_coorg_elem(myrank,npgeo_ori,coorg,numat,ngnod,nspec, & pointsdisp,plot_lowerleft_corner_only, & nelemabs,nelem_acoustic_surface, & num_fluid_solid_edges,num_fluid_poro_edges, & num_solid_poro_edges,nnodes_tangential_curve) !---- allocate arrays - if(ipass == 1) then allocate(shape2D(ngnod,NGLLX,NGLLZ)) allocate(dershape2D(NDIM,ngnod,NGLLX,NGLLZ)) allocate(shape2D_display(ngnod,pointsdisp,pointsdisp)) @@ -1234,25 +1191,19 @@ program specfem2D allocate(inv_tau_sigma_nu2_sent(N_SLS)) allocate(phi_nu1_sent(N_SLS)) allocate(phi_nu2_sent(N_SLS)) - endif ! !---- read the material properties ! call gmat01(density,porosity,tortuosity,anisotropy,permeability,poroelastcoef,numat,& - myrank,ipass,QKappa_attenuation,Qmu_attenuation,freq0,Q0,f0(1),ATTENUATION_PORO_FLUID_PART) + myrank,QKappa_attenuation,Qmu_attenuation,freq0,Q0,f0(1),ATTENUATION_PORO_FLUID_PART) ! !---- read spectral macrobloc data ! - if(ipass == 1) then - allocate(antecedent_list(nspec)) - allocate(perm(nspec)) - endif ! DK DK add support for using pml in mpi mode with external mesh allocate(region_CPML(nspec)) - call read_databases_mato(ipass,nspec,ngnod,kmato,knods, & - perm,antecedent_list,region_CPML) + call read_databases_mato(nspec,ngnod,kmato,knods,region_CPML) !! DK DK Dec 2011: add a small crack (discontinuity) in the medium manually if(ADD_A_SMALL_CRACK_IN_THE_MEDIUM) then @@ -1271,8 +1222,6 @@ program specfem2D if(ngnod /= 4) stop 'must currently have ngnod == 4 when adding a crack manually' - if(PERFORM_CUTHILL_MCKEE) stop 'must not have PERFORM_CUTHILL_MCKEE when adding a crack manually' - if(FAST_NUMBERING) stop 'must not have FAST_NUMBERING when adding a crack manually' !! DK DK modify arrays "knods" and "coorg" to introduce the crack manually by duplicating and splitting the nodes @@ -1366,7 +1315,6 @@ program specfem2D #endif ! allocate memory variables for attenuation - if(ipass == 1) then allocate(e1(NGLLX,NGLLZ,nspec_allocate,N_SLS)) allocate(e11(NGLLX,NGLLZ,nspec_allocate,N_SLS)) allocate(e13(NGLLX,NGLLZ,nspec_allocate,N_SLS)) @@ -1411,7 +1359,6 @@ program specfem2D e13_force_rk(:,:,:,:,:) = 0._CUSTOM_REAL allocate(Mu_nu1(NGLLX,NGLLZ,nspec)) allocate(Mu_nu2(NGLLX,NGLLZ,nspec)) - endif ! define the attenuation quality factors. ! they can be different for each element. @@ -1433,7 +1380,6 @@ program specfem2D enddo ! allocate memory variables for viscous attenuation (poroelastic media) - if(ipass == 1) then if(ATTENUATION_PORO_FLUID_PART) then allocate(rx_viscous(NGLLX,NGLLZ,nspec)) allocate(rz_viscous(NGLLX,NGLLZ,nspec)) @@ -1458,14 +1404,12 @@ program specfem2D allocate(viscox(NGLLX,NGLLZ,1)) allocate(viscoz(NGLLX,NGLLZ,1)) endif - endif ! !---- read interfaces data ! call read_databases_ninterface(ninterface,max_interface_size) if ( ninterface > 0 ) then - if(ipass == 1) then allocate(my_neighbours(ninterface)) allocate(my_nelmnts_neighbours(ninterface)) allocate(my_interfaces(4,max_interface_size,ninterface)) @@ -1478,13 +1422,9 @@ program specfem2D allocate(inum_interfaces_acoustic(ninterface)) allocate(inum_interfaces_elastic(ninterface)) allocate(inum_interfaces_poroelastic(ninterface)) - endif - call read_databases_interfaces(ipass,ninterface,nspec,max_interface_size, & - my_neighbours,my_nelmnts_neighbours,my_interfaces, & - perm,antecedent_list) + call read_databases_interfaces(ninterface,max_interface_size,my_neighbours,my_nelmnts_neighbours,my_interfaces) else - if(ipass == 1) then allocate(my_neighbours(1)) allocate(my_nelmnts_neighbours(1)) allocate(my_interfaces(1,1,1)) @@ -1497,7 +1437,6 @@ program specfem2D allocate(inum_interfaces_acoustic(1)) allocate(inum_interfaces_elastic(1)) allocate(inum_interfaces_poroelastic(1)) - endif endif @@ -1510,7 +1449,6 @@ program specfem2D anyabs = .true. endif - if(ipass == 1) then allocate(numabs(nelemabs)) allocate(codeabs(4,nelemabs)) allocate(typeabs(nelemabs)) @@ -1540,15 +1478,13 @@ program specfem2D allocate(ib_bottom(nelemabs)) allocate(ib_top(nelemabs)) - endif - ! !---- read absorbing boundary data ! - call read_databases_absorbing(myrank,ipass,nelemabs,nspec,anyabs, & + call read_databases_absorbing(myrank,nelemabs,nspec,anyabs, & ibegin_edge1,iend_edge1,ibegin_edge2,iend_edge2, & ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4, & - numabs,codeabs,typeabs,perm,antecedent_list, & + numabs,codeabs,typeabs, & nspec_left,nspec_right,nspec_bottom,nspec_top, & ib_right,ib_left,ib_bottom,ib_top,PML_BOUNDARY_CONDITIONS) @@ -1560,8 +1496,7 @@ program specfem2D if( anyabs ) then - ! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method - if(ipass == 1) then + ! files to save absorbed waves needed to reconstruct backward wavefield for adjoint method if(any_elastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE == 3).and. (.not. PML_BOUNDARY_CONDITIONS)) then allocate(b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP)) allocate(b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP)) @@ -1603,7 +1538,6 @@ program specfem2D allocate(b_absorb_acoustic_bottom(1,1,1)) allocate(b_absorb_acoustic_top(1,1,1)) endif - endif else @@ -1643,12 +1577,9 @@ program specfem2D any_acoustic_edges = .false. nelem_acoustic_surface = 1 endif - if( ipass == 1 ) then - allocate(acoustic_edges(4,nelem_acoustic_surface)) - allocate(acoustic_surface(5,nelem_acoustic_surface)) - endif - call read_databases_free_surf(ipass,nelem_acoustic_surface,nspec, & - acoustic_edges,perm,antecedent_list,any_acoustic_edges) + allocate(acoustic_edges(4,nelem_acoustic_surface)) + allocate(acoustic_surface(5,nelem_acoustic_surface)) + call read_databases_free_surf(nelem_acoustic_surface,acoustic_edges,any_acoustic_edges) ! resets nelem_acoustic_surface if( any_acoustic_edges .eqv. .false. ) nelem_acoustic_surface = 0 @@ -1656,7 +1587,7 @@ program specfem2D if(nelem_acoustic_surface > 0) then call construct_acoustic_surface ( nspec, ngnod, knods, nelem_acoustic_surface, & acoustic_edges, acoustic_surface) - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then write(IOUT,*) write(IOUT,*) 'Number of free surface elements: ',nelem_acoustic_surface endif @@ -1672,44 +1603,37 @@ program specfem2D any_fluid_solid_edges = .false. num_fluid_solid_edges = 1 endif - if(ipass == 1) then - allocate(fluid_solid_acoustic_ispec(num_fluid_solid_edges)) - allocate(fluid_solid_acoustic_iedge(num_fluid_solid_edges)) - allocate(fluid_solid_elastic_ispec(num_fluid_solid_edges)) - allocate(fluid_solid_elastic_iedge(num_fluid_solid_edges)) - endif + allocate(fluid_solid_acoustic_ispec(num_fluid_solid_edges)) + allocate(fluid_solid_acoustic_iedge(num_fluid_solid_edges)) + allocate(fluid_solid_elastic_ispec(num_fluid_solid_edges)) + allocate(fluid_solid_elastic_iedge(num_fluid_solid_edges)) if( num_fluid_poro_edges > 0 ) then any_fluid_poro_edges = .true. else any_fluid_poro_edges = .false. num_fluid_poro_edges = 1 endif - if(ipass == 1) then - allocate(fluid_poro_acoustic_ispec(num_fluid_poro_edges)) - allocate(fluid_poro_acoustic_iedge(num_fluid_poro_edges)) - allocate(fluid_poro_poroelastic_ispec(num_fluid_poro_edges)) - allocate(fluid_poro_poroelastic_iedge(num_fluid_poro_edges)) - endif + allocate(fluid_poro_acoustic_ispec(num_fluid_poro_edges)) + allocate(fluid_poro_acoustic_iedge(num_fluid_poro_edges)) + allocate(fluid_poro_poroelastic_ispec(num_fluid_poro_edges)) + allocate(fluid_poro_poroelastic_iedge(num_fluid_poro_edges)) if ( num_solid_poro_edges > 0 ) then any_solid_poro_edges = .true. else any_solid_poro_edges = .false. num_solid_poro_edges = 1 endif - if(ipass == 1) then - allocate(solid_poro_elastic_ispec(num_solid_poro_edges)) - allocate(solid_poro_elastic_iedge(num_solid_poro_edges)) - allocate(solid_poro_poroelastic_ispec(num_solid_poro_edges)) - allocate(solid_poro_poroelastic_iedge(num_solid_poro_edges)) - endif + allocate(solid_poro_elastic_ispec(num_solid_poro_edges)) + allocate(solid_poro_elastic_iedge(num_solid_poro_edges)) + allocate(solid_poro_poroelastic_ispec(num_solid_poro_edges)) + allocate(solid_poro_poroelastic_iedge(num_solid_poro_edges)) - call read_databases_coupled(ipass,nspec,num_fluid_solid_edges,any_fluid_solid_edges, & + call read_databases_coupled(num_fluid_solid_edges,any_fluid_solid_edges, & fluid_solid_acoustic_ispec,fluid_solid_elastic_ispec, & num_fluid_poro_edges,any_fluid_poro_edges, & fluid_poro_acoustic_ispec,fluid_poro_poroelastic_ispec, & num_solid_poro_edges,any_solid_poro_edges, & - solid_poro_elastic_ispec,solid_poro_poroelastic_ispec, & - perm,antecedent_list) + solid_poro_elastic_ispec,solid_poro_poroelastic_ispec) ! resets counters if( any_fluid_solid_edges .eqv. .false. ) num_fluid_solid_edges = 0 @@ -1727,10 +1651,8 @@ program specfem2D any_tangential_curve = .false. nnodes_tangential_curve = 1 endif - if (ipass == 1) then - allocate(nodes_tangential_curve(2,nnodes_tangential_curve)) - allocate(dist_tangential_detection_curve(nnodes_tangential_curve)) - endif + allocate(nodes_tangential_curve(2,nnodes_tangential_curve)) + allocate(dist_tangential_detection_curve(nnodes_tangential_curve)) call read_databases_final(nnodes_tangential_curve,nodes_tangential_curve, & force_normal_to_surface,rec_normal_to_surface, & any_tangential_curve) @@ -1756,9 +1678,9 @@ program specfem2D ! "slow and clean" or "quick and dirty" version if(FAST_NUMBERING) then - call createnum_fast(knods,ibool,shape2D,coorg,nglob,npgeo,nspec,ngnod,myrank,ipass) + call createnum_fast(knods,ibool,shape2D,coorg,nglob,npgeo,nspec,ngnod,myrank) else - call createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank,ipass) + call createnum_slow(knods,ibool,nglob,nspec,ngnod,myrank) endif #ifdef USE_MPI @@ -1770,7 +1692,7 @@ program specfem2D nspec_total = nspec nglob_total = nglob #endif - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then write(IOUT,*) write(IOUT,*) 'Total number of elements: ',nspec_total write(IOUT,*) 'decomposed as follows:' @@ -1808,11 +1730,6 @@ program specfem2D endif endif -! create a new indirect addressing array to reduce cache misses in memory access in the solver - if(ipass == 2) then - - deallocate(perm) - !! DK DK for periodic conditions: detect common points between left and right edges if(ADD_PERIODIC_CONDITIONS) then @@ -1824,9 +1741,6 @@ program specfem2D if(any_poroelastic .or. any_acoustic) & stop 'periodic conditions currently implemented for purely elastic models only' - if(ACTUALLY_IMPLEMENT_PERM_OUT .or. ACTUALLY_IMPLEMENT_PERM_INN .or. ACTUALLY_IMPLEMENT_PERM_WHOLE) & - stop 'currently, all permutations should be off for periodic conditions' - print * open(unit=123,file='DATA/Database00000_left_edge_only',status='old') read(123,*) NSPEC_PERIO @@ -1957,13 +1871,7 @@ program specfem2D !! DK DK end of periodic conditions: detect common points between left and right edges ! reduces cache misses - call get_global(nspec_outer,nspec,nglob,ibool) - - else if(ipass /= 1) then - - stop 'incorrect pass number for reduction of cache misses' - - endif ! ipass + call get_global(nspec,nglob,ibool) !---- compute shape functions and their derivatives for regular interpolated display grid do j = 1,pointsdisp @@ -1992,7 +1900,7 @@ program specfem2D enddo close(IIN) - if (myrank == 0 .and. ipass == 1) then + if (myrank == 0) then write(IOUT,*) write(IOUT,*) 'Total number of receivers = ',nrec write(IOUT,*) @@ -2001,8 +1909,6 @@ program specfem2D if(nrec < 1) call exit_MPI('need at least one receiver') ! receiver information - if(ipass == 1) then - allocate(ispec_selected_rec(nrec)) allocate(st_xval(nrec)) allocate(st_zval(nrec)) @@ -2073,8 +1979,6 @@ program specfem2D allocate(c25ext(1,1,1)) endif - endif - ! !---- set the coordinates of the points of the global grid ! @@ -2174,7 +2078,7 @@ program specfem2D ! !--- save the grid of points in a file ! - if(output_grid_ASCII .and. myrank == 0 .and. ipass == 1) then + if(output_grid_ASCII .and. myrank == 0) then write(IOUT,*) write(IOUT,*) 'Saving the grid in an ASCII text file...' write(IOUT,*) @@ -2189,12 +2093,11 @@ program specfem2D ! !----- plot the GLL mesh in a Gnuplot file ! - if(output_grid_Gnuplot .and. myrank == 0 .and. ipass == 1) & + if(output_grid_Gnuplot .and. myrank == 0) & call plotgll(knods,ibool,coorg,coord,nglob,npgeo,ngnod,nspec) -! if (assign_external_model .and. ipass == 1) then if (assign_external_model) then - if(myrank == 0 .and. ipass == 1) write(IOUT,*) 'Assigning an external velocity and density model...' + if(myrank == 0) write(IOUT,*) 'Assigning an external velocity and density model...' call read_external_model(any_acoustic,any_elastic,any_poroelastic, & elastic,poroelastic,anisotropic,nspec,nglob,N_SLS,ibool, & f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, & @@ -2256,17 +2159,17 @@ program specfem2D !---- define actual location of source and receivers - call setup_sources_receivers(NSOURCES,initialfield,source_type,& + call setup_sources_receivers(NSOURCES,initialfield,source_type, & coord,ibool,nglob,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, & x_source,z_source,ispec_selected_source,ispec_selected_rec, & - is_proc_source,nb_proc_source,ipass,& - sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,& + is_proc_source,nb_proc_source, & + sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo, & nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod, & nrec,nrecloc,recloc,which_proc_receiver,st_xval,st_zval, & xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver,iglob_source) ! compute source array for adjoint source - if(ipass == 1) nadj_rec_local = 0 + nadj_rec_local = 0 if(SIMULATION_TYPE == 3) then ! adjoint calculation do irec = 1,nrec @@ -2277,10 +2180,10 @@ program specfem2D nadj_rec_local = nadj_rec_local + 1 endif enddo - if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ)) - if (nadj_rec_local > 0 .and. ipass == 1) then + allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ)) + if (nadj_rec_local > 0) then allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ)) - else if (ipass == 1) then + else allocate(adj_sourcearrays(1,1,1,1,1)) endif @@ -2332,11 +2235,10 @@ program specfem2D close(113) deallocate(adj_src_s) endif - else if (ipass == 1) then + else allocate(adj_sourcearrays(1,1,1,1,1)) endif - if (ipass == 1) then if (nrecloc > 0) then allocate(anglerec_irec(nrecloc)) allocate(cosrot_irec(nrecloc)) @@ -2355,12 +2257,10 @@ program specfem2D anglerec_irec(:) = anglerec * pi / 180.d0 cosrot_irec(:) = cos(anglerec_irec(:)) sinrot_irec(:) = sin(anglerec_irec(:)) - endif ! !--- tangential computation ! - if (ipass == NUMBER_OF_PASSES) then ! for receivers if (rec_normal_to_surface) then @@ -2549,18 +2449,14 @@ program specfem2D endif ! force_normal_to_surface - endif ! ipass - ! !--- ! ! allocate seismogram arrays - if(ipass == 1) then - allocate(sisux(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc)) - allocate(sisuz(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc)) - allocate(siscurl(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc)) - endif + allocate(sisux(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc)) + allocate(sisuz(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc)) + allocate(siscurl(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,nrecloc)) ! check if acoustic receiver is exactly on the free surface because pressure is zero there do ispec_acoustic_surface = 1,nelem_acoustic_surface @@ -2620,8 +2516,6 @@ program specfem2D enddo ! displacement, velocity, acceleration and inverse of the mass matrix for elastic elements - if(ipass == 1) then - if(any_elastic) then nglob_elastic = nglob else @@ -3233,7 +3127,6 @@ program specfem2D allocate(alpha_x_store(1,1,1)) allocate(alpha_z_store(1,1,1)) endif ! PML_BOUNDARY_CONDITIONS - endif ! ipass == 1 ! !---- build the global mass matrix @@ -3261,7 +3154,7 @@ program specfem2D if ( nproc > 1 ) then ! preparing for MPI communications - if(ipass == 1) allocate(mask_ispec_inner_outer(nspec)) + allocate(mask_ispec_inner_outer(nspec)) mask_ispec_inner_outer(:) = .false. call get_MPI(nspec,ibool,knods,ngnod,nglob,elastic,poroelastic, & @@ -3275,19 +3168,15 @@ program specfem2D inum_interfaces_poroelastic, & ninterface_acoustic, ninterface_elastic, ninterface_poroelastic, & mask_ispec_inner_outer, & - myrank,ipass,coord) - + myrank,coord) nspec_outer = count(mask_ispec_inner_outer) nspec_inner = nspec - nspec_outer - if(ipass == 1) then - allocate(ispec_outer_to_glob(nspec_outer)) - allocate(ispec_inner_to_glob(nspec_inner)) - endif + allocate(ispec_outer_to_glob(nspec_outer)) + allocate(ispec_inner_to_glob(nspec_inner)) ! building of corresponding arrays between inner/outer elements and their global number - if(ipass == 1) then num_ispec_outer = 0 num_ispec_inner = 0 do ispec = 1, nspec @@ -3299,13 +3188,11 @@ program specfem2D ispec_inner_to_glob(num_ispec_inner) = ispec endif enddo - endif ! buffers for MPI communications max_ibool_interfaces_size_ac = maxval(nibool_interfaces_acoustic(:)) max_ibool_interfaces_size_el = 3*maxval(nibool_interfaces_elastic(:)) max_ibool_interfaces_size_po = NDIM*maxval(nibool_interfaces_poroelastic(:)) - if(ipass == 1) then allocate(tab_requests_send_recv_acoustic(ninterface_acoustic*2)) allocate(buffer_send_faces_vector_ac(max_ibool_interfaces_size_ac,ninterface_acoustic)) allocate(buffer_recv_faces_vector_ac(max_ibool_interfaces_size_ac,ninterface_acoustic)) @@ -3317,7 +3204,6 @@ program specfem2D allocate(buffer_recv_faces_vector_pos(max_ibool_interfaces_size_po,ninterface_poroelastic)) allocate(buffer_send_faces_vector_pow(max_ibool_interfaces_size_po,ninterface_poroelastic)) allocate(buffer_recv_faces_vector_pow(max_ibool_interfaces_size_po,ninterface_poroelastic)) - endif ! assembling the mass matrix call assemble_MPI_scalar(rmass_inverse_acoustic,nglob_acoustic, & @@ -3338,52 +3224,44 @@ program specfem2D num_ispec_outer = 0 num_ispec_inner = 0 - if(ipass == 1) allocate(mask_ispec_inner_outer(1)) + allocate(mask_ispec_inner_outer(1)) nspec_outer = 0 nspec_inner = nspec - if(ipass == 1) allocate(ispec_inner_to_glob(nspec_inner)) + allocate(ispec_inner_to_glob(nspec_inner)) do ispec = 1, nspec ispec_inner_to_glob(ispec) = ispec enddo - endif ! end of test on wether there is more than one process (nproc > 1) + endif ! end of test on whether there is more than one process (nproc > 1) #else num_ispec_outer = 0 num_ispec_inner = 0 - if(ipass == 1) allocate(mask_ispec_inner_outer(1)) + allocate(mask_ispec_inner_outer(1)) nspec_outer = 0 nspec_inner = nspec - if(ipass == 1) then - allocate(ispec_outer_to_glob(1)) - allocate(ispec_inner_to_glob(nspec_inner)) - endif + allocate(ispec_outer_to_glob(1)) + allocate(ispec_inner_to_glob(nspec_inner)) do ispec = 1, nspec ispec_inner_to_glob(ispec) = ispec enddo #endif - if(ipass == 1) then - - ! allocate(antecedent_list(nspec)) - ! loop over spectral elements do ispec_outer = 1,nspec_outer ! get global numbering for inner or outer elements ispec = ispec_outer_to_glob(ispec_outer) - antecedent_list(ispec) = ispec_outer enddo ! loop over spectral elements do ispec_inner = 1,nspec_inner ! get global numbering for inner or outer elements ispec = ispec_inner_to_glob(ispec_inner) - antecedent_list(ispec) = nspec_outer + ispec_inner enddo allocate(ibool_outer(NGLLX,NGLLZ,nspec_outer)) @@ -3412,85 +3290,15 @@ program specfem2D ! reduces cache misses for inner elements call get_global_indirect_addressing(nspec_inner,nglob,ibool_inner) - ! the total number of points without multiples in this region is now known nglob_inner = maxval(ibool_inner) - !allocate(perm(nspec)) - - ! use identity permutation by default - do ispec = 1,nspec - perm(ispec) = ispec - enddo - - if(ACTUALLY_IMPLEMENT_PERM_WHOLE) then - - allocate(check_perm(nspec)) - call get_perm_cuthill_mckee(ibool,perm,LIMIT_MULTI_CUTHILL,nspec,nglob,PERFORM_CUTHILL_MCKEE) - ! check that the permutation obtained is bijective - check_perm(:) = -1 - do ispec = 1,nspec - check_perm(perm(ispec)) = ispec - enddo - if(minval(check_perm) /= 1) stop 'minval check_perm is incorrect for whole' - if(maxval(check_perm) /= nspec) stop 'maxval check_perm is incorrect for whole' - deallocate(check_perm) - else - - if(ACTUALLY_IMPLEMENT_PERM_OUT) then - allocate(check_perm(nspec_outer)) - call get_perm_cuthill_mckee(ibool_outer,perm(1:nspec_outer),LIMIT_MULTI_CUTHILL,nspec_outer,nglob_outer,PERFORM_CUTHILL_MCKEE) - ! check that the permutation obtained is bijective - check_perm(:) = -1 - do ispec = 1,nspec_outer - check_perm(perm(ispec)) = ispec - enddo - if(minval(check_perm) /= 1) stop 'minval check_perm is incorrect for outer' - if(maxval(check_perm) /= nspec_outer) stop 'maxval check_perm is incorrect for outer' - deallocate(check_perm) - deallocate(ibool_outer) - endif - - if(ACTUALLY_IMPLEMENT_PERM_INN) then - allocate(check_perm(nspec_inner)) - call get_perm_cuthill_mckee(ibool_inner,perm(nspec_outer+1:nspec),LIMIT_MULTI_CUTHILL,nspec_inner,nglob_inner, & - PERFORM_CUTHILL_MCKEE) - ! check that the permutation obtained is bijective - check_perm(:) = -1 - do ispec = 1,nspec_inner - check_perm(perm(nspec_outer+ispec)) = ispec - enddo - if(minval(check_perm) /= 1) stop 'minval check_perm is incorrect for inner' - if(maxval(check_perm) /= nspec_inner) stop 'maxval check_perm is incorrect for inner' - deallocate(check_perm) - ! add the right offset - perm(nspec_outer+1:nspec) = perm(nspec_outer+1:nspec) + nspec_outer - deallocate(ibool_inner) - endif - - endif - - endif - - enddo ! end of further reduction of cache misses inner/outer in two passes - -!============================================ -! -! end inner/outer passes -! -!============================================ - -!--- -!--- end of section performed in two passes -!--- - call invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,& rmass_inverse_elastic_one,rmass_inverse_elastic_three,nglob_elastic, & rmass_inverse_acoustic,nglob_acoustic, & rmass_s_inverse_poroelastic, & rmass_w_inverse_poroelastic,nglob_poroelastic) - ! check the mesh, stability and number of points per wavelength if(DISPLAY_SUBSET_OPTION == 1) then UPPER_LIMIT_DISPLAY = nspec @@ -4889,11 +4697,6 @@ program specfem2D ! to dump the wave field this_is_the_first_time_we_dump = .true. -#ifdef USE_MPI -! add a barrier if we generate traces of the run for analysis with "ParaVer" - if(GENERATE_PARAVER_TRACES) call MPI_BARRIER(MPI_COMM_WORLD,ier) -#endif - ! !---- s t a r t t i m e i t e r a t i o n s ! @@ -7759,10 +7562,6 @@ program specfem2D any_acoustic_glob,any_acoustic,potential_acoustic, & timestamp_seconds_start) -#ifdef USE_MPI -! add a barrier if we generate traces of the run for analysis with "ParaVer" - if(GENERATE_PARAVER_TRACES) call MPI_BARRIER(MPI_COMM_WORLD,ier) -#endif endif !---- loop on all the receivers to compute and store the seismograms @@ -8866,7 +8665,6 @@ program specfem2D !---- save temporary or final seismograms ! suppress seismograms if we generate traces of the run for analysis with "ParaVer", because time consuming if(mod(it,NSTEP_BETWEEN_OUTPUT_SEISMOS) == 0 .or. it == NSTEP) then - if(.not. GENERATE_PARAVER_TRACES) & call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, & nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, & NSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current,p_sv, &