From deae98344e62308424b29109dc41ec88bb46b3d5 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Tue, 2 Jul 2019 15:27:38 +0200 Subject: [PATCH] code cleaning (removing unused parameters; using modules instead of include file); fix in reading poroelastic file lines; increasing precision for irregular element checks to avoid underflows --- src/decompose_mesh/fault_scotch.f90 | 2 + src/decompose_mesh/lts_helper.F90 | 2 +- src/decompose_mesh/lts_setup_elements.F90 | 3 + src/decompose_mesh/partition_scotch.F90 | 406 ++++--- src/decompose_mesh/read_mesh_files.F90 | 9 +- src/decompose_mesh/rules.mk | 1 + src/generate_databases/calc_jacobian.f90 | 17 +- .../create_mass_matrices.f90 | 10 +- .../create_regions_mesh.f90 | 2 +- .../fault_generate_databases.f90 | 12 +- src/generate_databases/get_model.F90 | 3 +- src/generate_databases/model_coupled.f90 | 2 +- src/generate_databases/model_sep.f90 | 12 +- .../pml_set_local_dampingcoeff.f90 | 4 +- src/generate_databases/save_arrays_solver.F90 | 4 +- src/generate_databases/setup_color_perm.f90 | 2 +- .../input_output/IO_model_mod.f90 | 4 +- .../input_output/mesh_tools_mod.f90 | 3 +- src/meshfem3D/check_mesh_quality.f90 | 3 +- src/meshfem3D/compute_parameters.f90 | 3 +- src/meshfem3D/create_interfaces_mesh.f90 | 2 +- src/meshfem3D/create_meshfem_mesh.f90 | 1064 +++++++++-------- src/meshfem3D/create_visual_files.f90 | 3 +- src/meshfem3D/define_subregions_heuristic.f90 | 3 +- src/meshfem3D/define_superbrick.f90 | 9 +- src/meshfem3D/determine_cavity.f90 | 10 +- src/meshfem3D/earth_chunk.f90 | 4 +- src/meshfem3D/get_MPI_cutplanes_eta.f90 | 3 +- src/meshfem3D/get_MPI_cutplanes_xi.f90 | 3 +- src/meshfem3D/get_flags_boundaries.f90 | 3 +- src/meshfem3D/meshfem3D_adios_stubs.f90 | 4 +- src/meshfem3D/rules.mk | 20 +- src/meshfem3D/save_databases.F90 | 18 +- src/meshfem3D/save_databases_adios.F90 | 3 +- src/meshfem3D/store_boundaries.f90 | 9 +- src/meshfem3D/store_coords.f90 | 3 +- src/shared/get_attenuation_model.f90 | 2 +- src/shared/gll_library.f90 | 16 +- src/shared/parallel.f90 | 2 - src/shared/read_parameter_file.F90 | 2 +- src/specfem3D/compute_stacey_viscoelastic.F90 | 2 +- src/specfem3D/get_cmt.f90 | 2 +- src/specfem3D/get_force.f90 | 2 +- src/specfem3D/make_gravity.f90 | 196 +-- .../pml_compute_accel_contribution.f90 | 23 +- .../pml_compute_memory_variables.f90 | 23 +- src/specfem3D/prepare_gpu.f90 | 310 ++--- src/specfem3D/prepare_timerun.F90 | 2 +- src/specfem3D/read_stations.f90 | 2 +- .../smooth_sem.F90 | 4 +- src/tomography/read_parameters_tomo.f90 | 17 +- 51 files changed, 1164 insertions(+), 1106 deletions(-) diff --git a/src/decompose_mesh/fault_scotch.f90 b/src/decompose_mesh/fault_scotch.f90 index 70ffaa8b7..914c30db2 100644 --- a/src/decompose_mesh/fault_scotch.f90 +++ b/src/decompose_mesh/fault_scotch.f90 @@ -71,7 +71,9 @@ subroutine read_fault_files(localpath_name) read(101,*) nbfaults else nbfaults = 0 + print * print *, 'Par_file_faults not found: assuming that there are no faults' + print * endif close(101) diff --git a/src/decompose_mesh/lts_helper.F90 b/src/decompose_mesh/lts_helper.F90 index da58dfcaa..b75e4edf2 100644 --- a/src/decompose_mesh/lts_helper.F90 +++ b/src/decompose_mesh/lts_helper.F90 @@ -100,7 +100,7 @@ end subroutine lts_add_overlap_elements subroutine lts_save_p_level_partitions() - use constants, only: IOVTK,MAX_STRING_LEN + use constants, only: MAX_STRING_LEN use decompose_mesh_par, only: part, ispec_p_refine, nspec, nnodes, nodes_coords, elmnts, NGNOD, LOCAL_PATH diff --git a/src/decompose_mesh/lts_setup_elements.F90 b/src/decompose_mesh/lts_setup_elements.F90 index a721decc8..7b14d09d5 100644 --- a/src/decompose_mesh/lts_setup_elements.F90 +++ b/src/decompose_mesh/lts_setup_elements.F90 @@ -89,9 +89,12 @@ subroutine lts_setup_elements() ! user output if (LTS_MODE) then + print * print *,'local time stepping: turned ON' else + print * print *,'local time stepping: turned OFF' + print * ! nothing to do return endif diff --git a/src/decompose_mesh/partition_scotch.F90 b/src/decompose_mesh/partition_scotch.F90 index d6a654503..8d0a7028f 100644 --- a/src/decompose_mesh/partition_scotch.F90 +++ b/src/decompose_mesh/partition_scotch.F90 @@ -55,6 +55,25 @@ ! To select SCOTCH partioning, in Par_file choose: ! PARTITIONING_TYPE = 1 +module scotch_par + + implicit none + +#if defined(USE_SCOTCH) + include 'scotchf.h' + + ! SCOTCH + double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph + double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat +!!!!!! character(len=*), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}' + +#endif + +end module scotch_par + +! +!---------------------------------------------------------------------------------------------- +! subroutine partition_scotch() @@ -65,25 +84,13 @@ subroutine partition_scotch() use decompose_mesh_par, only: NSPEC, nparts, part, nb_edges, elmnts_load, xadj, adjncy, & LTS_MODE - implicit none + use scotch_par - include 'scotchf.h' + implicit none ! local parameters integer :: ier - ! scotch p-level partitioning - integer,dimension(:),allocatable :: part_tmp,elmnts_load_tmp - integer,dimension(:),allocatable :: ispec_global,ispec_local - integer,dimension(:),allocatable :: xadj_tmp,adjncy_tmp - integer :: nb_edges_tmp - integer :: i,j,k,ispec,iproc - - ! SCOTCH - double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph - double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat -!!!!!! character(len=*), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}' - ! checks if partitioning needed if (nparts == 1) return @@ -149,6 +156,7 @@ subroutine partition_scotch() if (LTS_MODE .and. BALANCE_P_LEVELS) then ! LTS mode partitioning for each p-level call lts_partition_scotch() + else ! over-all scotch partitioning call scotchfgraphbuild (scotchgraph (1), 0, nspec, & @@ -182,208 +190,226 @@ subroutine partition_scotch() stop 'Scotch ERROR : MAIN : Cannot destroy strategy' endif - contains +#else - subroutine lts_partition_scotch() + ! compiled without support + print *,'Scotch Error: Please re-compile with -DUSE_SCOTCH to use SCOTCH partitioning.' + stop 'No SCOTCH support' - ! note: SCOTCH is unable to handle constrained partitioning as Patoh. - ! the goal would be to have about the same number of p-elements in each partition (due to MPI scheme). - ! as an ugly work-around, we partition for each p-level again, using newly assigned element weights +#endif - use constants, only: BALANCE_P_LEVELS_PARTIAL,P_LEVEL_PARTIAL_SUBSET_MINIMUM + end subroutine partition_scotch - use decompose_mesh_par, only: p_level, num_p_level, ispec_p_refine, num_ispec_level +! +!---------------------------------------------------------------------------------------------- +! - implicit none +#if defined(USE_SCOTCH) - integer :: ilevel,nspec_p - integer :: nparts_partial - integer :: p,inum + subroutine lts_partition_scotch() - ! allocates temporary arrays - allocate(part_tmp(1:nspec), & - ispec_global(1:nspec), & - ispec_local(1:nspec), & - elmnts_load_tmp(1:nspec), & - xadj_tmp(1:nspec+1), & - adjncy_tmp(1:nb_edges+1), & - stat=ier) - if (ier /= 0) stop 'Error allocating array part_tmp,...' + ! note: SCOTCH is unable to handle constrained partitioning as Patoh. + ! the goal would be to have about the same number of p-elements in each partition (due to MPI scheme). + ! as an ugly work-around, we partition for each p-level again, using newly assigned element weights - ! partitions each p-level separately - do ilevel = 1,num_p_level - p = p_level(ilevel) - nspec_p = num_ispec_level(ilevel) + use decompose_mesh_par, only: NSPEC, nparts, part, nb_edges, elmnts_load, xadj, adjncy + ! p-levels + use constants, only: BALANCE_P_LEVELS_PARTIAL,P_LEVEL_PARTIAL_SUBSET_MINIMUM + use decompose_mesh_par, only: p_level, num_p_level, ispec_p_refine, num_ispec_level - ! element lookup tables - ispec_global(:) = 0 - ispec_local(:) = 0 - inum = 0 - do ispec = 1,nspec - if (ispec_p_refine(ispec) == p) then - inum = inum + 1 - if (inum > nspec_p) stop 'Error index out of bounds of nspec_p' - ispec_global(inum) = ispec - ispec_local(ispec) = inum + use scotch_par + + implicit none + + ! local parameters + integer :: ier + + ! scotch p-level partitioning + integer,dimension(:),allocatable :: part_tmp,elmnts_load_tmp + integer,dimension(:),allocatable :: ispec_global,ispec_local + integer,dimension(:),allocatable :: xadj_tmp,adjncy_tmp + integer :: nb_edges_tmp + integer :: i,j,k,ispec,iproc + + integer :: ilevel,nspec_p + integer :: nparts_partial + integer :: p,inum + + ! allocates temporary arrays + allocate(part_tmp(1:nspec), & + ispec_global(1:nspec), & + ispec_local(1:nspec), & + elmnts_load_tmp(1:nspec), & + xadj_tmp(1:nspec+1), & + adjncy_tmp(1:nb_edges+1), & + stat=ier) + if (ier /= 0) stop 'Error allocating array part_tmp,...' + + ! partitions each p-level separately + do ilevel = 1,num_p_level + p = p_level(ilevel) + nspec_p = num_ispec_level(ilevel) + + ! element lookup tables + ispec_global(:) = 0 + ispec_local(:) = 0 + inum = 0 + do ispec = 1,nspec + if (ispec_p_refine(ispec) == p) then + inum = inum + 1 + if (inum > nspec_p) stop 'Error index out of bounds of nspec_p' + ispec_global(inum) = ispec + ispec_local(ispec) = inum + endif + enddo + if (inum /= nspec_p) stop 'Error number of p-vertices not equals to number of p-elements' + + ! builds up local arrays only for p-elements + xadj_tmp(:) = 0 + adjncy_tmp(:) = -1 + elmnts_load_tmp(:) = 0 + nb_edges_tmp = 1 + do inum = 1,nspec_p + ispec = ispec_global(inum) + + ! element load + elmnts_load_tmp(inum) = elmnts_load(ispec) + + ! number of neighbors + ! note: xadj() values start from 0; adjncy() values start from 0 to nspec-1 + ! we shift these by +1 since arrays in this subroutine are defined between 1 to .. + ! + ! xadj(i) -> adjncy( . ) : values in xadj point to the first index in adjncy() for vertex i + ! adjncy() holds all indices of neighbor vertices(=elements) + ! note: for xadj_tmp and adjncy_tmp we start indexing at 1 + xadj_tmp(inum) = nb_edges_tmp + do j = xadj(ispec),xadj(ispec+1)-1 + ! gets neighbor index + k = adjncy(j+1) + 1 + ! checks bounds + if (k < 1 .or. k > nspec) stop 'Error adjncy index' + ! adds vertex if it belongs to p-elements + if (ispec_p_refine(k) == p) then + ! we store the local index between 1 and nspec_p + i = ispec_local( k ) + adjncy_tmp(nb_edges_tmp) = i + nb_edges_tmp = nb_edges_tmp + 1 endif enddo - if (inum /= nspec_p) stop 'Error number of p-vertices not equals to number of p-elements' - - ! builds up local arrays only for p-elements - xadj_tmp(:) = 0 - adjncy_tmp(:) = -1 - elmnts_load_tmp(:) = 0 - nb_edges_tmp = 1 - do inum = 1,nspec_p - ispec = ispec_global(inum) - - ! element load - elmnts_load_tmp(inum) = elmnts_load(ispec) - - ! number of neighbors - ! note: xadj() values start from 0; adjncy() values start from 0 to nspec-1 - ! we shift these by +1 since arrays in this subroutine are defined between 1 to .. - ! - ! xadj(i) -> adjncy( . ) : values in xadj point to the first index in adjncy() for vertex i - ! adjncy() holds all indices of neighbor vertices(=elements) - ! note: for xadj_tmp and adjncy_tmp we start indexing at 1 - xadj_tmp(inum) = nb_edges_tmp - do j = xadj(ispec),xadj(ispec+1)-1 - ! gets neighbor index - k = adjncy(j+1) + 1 - ! checks bounds - if (k < 1 .or. k > nspec) stop 'Error adjncy index' - ! adds vertex if it belongs to p-elements - if (ispec_p_refine(k) == p) then - ! we store the local index between 1 and nspec_p - i = ispec_local( k ) - adjncy_tmp(nb_edges_tmp) = i - nb_edges_tmp = nb_edges_tmp + 1 - endif - enddo - enddo - ! last entry for xadj for a contiguous range of indices ("compact edge array") - xadj_tmp(nspec_p+1) = nb_edges_tmp - nb_edges_tmp = nb_edges_tmp - 1 - - ! debug - !print *,'xadj:',xadj(ispec_global(1):ispec_global(1)+1),xadj(nspec-1:nspec+1) - !print *,'xadj:',xadj(nspec+1) - !print *,'adjncy:' - !do j = xadj(ispec_global(1)),xadj(ispec_global(1)+1) - ! print *,j-xadj(ispec_global(1))+1,j,adjncy(j+1)+1,ispec_local(adjncy(j+1)+1) - !enddo - !print * - !print *,'temporary:',ispec_global(1) - !print *,'xadj:',xadj_tmp(1:2),xadj_tmp(nspec_p-1:nspec_p+2) - !print *,'xadj:',xadj_tmp(nspec_p+1) - !print *,'adjncy:' - !do j = xadj_tmp(1),xadj_tmp(2) - ! print *,j-xadj_tmp(1)+1,j,adjncy_tmp(j),ispec_global(adjncy_tmp(j)) - !enddo - - ! checks ranges - if (minval(adjncy_tmp(1:nb_edges_tmp)) < 1 .or. maxval(adjncy_tmp(1:nb_edges_tmp)) > nspec_p) then - print *,'Error adjncy bounds invalid', minval(adjncy_tmp(1:nb_edges_tmp)),maxval(adjncy_tmp(1:nb_edges_tmp)) - stop 'Error adjncy bounds invalid' - endif + enddo + ! last entry for xadj for a contiguous range of indices ("compact edge array") + xadj_tmp(nspec_p+1) = nb_edges_tmp + nb_edges_tmp = nb_edges_tmp - 1 + + ! debug + !print *,'xadj:',xadj(ispec_global(1):ispec_global(1)+1),xadj(nspec-1:nspec+1) + !print *,'xadj:',xadj(nspec+1) + !print *,'adjncy:' + !do j = xadj(ispec_global(1)),xadj(ispec_global(1)+1) + ! print *,j-xadj(ispec_global(1))+1,j,adjncy(j+1)+1,ispec_local(adjncy(j+1)+1) + !enddo + !print * + !print *,'temporary:',ispec_global(1) + !print *,'xadj:',xadj_tmp(1:2),xadj_tmp(nspec_p-1:nspec_p+2) + !print *,'xadj:',xadj_tmp(nspec_p+1) + !print *,'adjncy:' + !do j = xadj_tmp(1),xadj_tmp(2) + ! print *,j-xadj_tmp(1)+1,j,adjncy_tmp(j),ispec_global(adjncy_tmp(j)) + !enddo + + ! checks ranges + if (minval(adjncy_tmp(1:nb_edges_tmp)) < 1 .or. maxval(adjncy_tmp(1:nb_edges_tmp)) > nspec_p) then + print *,'Error adjncy bounds invalid', minval(adjncy_tmp(1:nb_edges_tmp)),maxval(adjncy_tmp(1:nb_edges_tmp)) + stop 'Error adjncy bounds invalid' + endif - ! sets up p-element partitioning - part_tmp(:) = -1 + ! sets up p-element partitioning + part_tmp(:) = -1 + + ! user output + print *,' p-level partitioning: p-level =',ilevel,'p =',p,'elements =',nspec_p !,'edges =',nb_edges_tmp + + ! scotch partitioning + ! arguments: #(1) graph_structure #(2)baseval (either 0/1) #(3)vertnbr (number_of_vertices) + ! #(4)verttab (adjacency_index_array) #(5)vendtab (adjacency_end_index_array (optional)) + ! #(6)velotab (vertex_load_array (optional)) #(7)vlbltab (vertex_label_array) + ! #(7)edgenbr (number_of_arcs) #(8)edgetab (adjacency_array) + ! #(9)edlotab (arc_load_array (optional)) #(10)ierror + ! + ! Since, in Fortran, there is no null reference, passing the scotchfgraphbuild routine + ! a reference equal to verttab in the velotab or vlbltab fields makes them be considered as missing arrays. + ! The same holds for edlotab when it is passed a reference equal to edgetab. + ! Setting vendtab to refer to one cell after verttab yields the same result, + ! as it is the exact semantics of a compact vertex array. + + call scotchfgraphbuild (scotchgraph(1), 1, nspec_p, & + xadj_tmp(1), xadj_tmp(1), & + elmnts_load_tmp(1), xadj_tmp(1), & + nb_edges_tmp, adjncy_tmp(1), & + adjncy_tmp(1), ier) + if (ier /= 0) stop 'ERROR : MAIN : Cannot build graph' + + ! checks scotch graph + call scotchfgraphcheck (scotchgraph(1), ier) + if (ier /= 0) stop 'Scotch ERROR : MAIN : Invalid check' - ! user output - print *,' p-level partitioning: p-level =',ilevel,'p =',p,'elements =',nspec_p !,'edges =',nb_edges_tmp - - ! scotch partitioning - ! arguments: #(1) graph_structure #(2)baseval (either 0/1) #(3)vertnbr (number_of_vertices) - ! #(4)verttab (adjacency_index_array) #(5)vendtab (adjacency_end_index_array (optional)) - ! #(6)velotab (vertex_load_array (optional)) #(7)vlbltab (vertex_label_array) - ! #(7)edgenbr (number_of_arcs) #(8)edgetab (adjacency_array) - ! #(9)edlotab (arc_load_array (optional)) #(10)ierror - ! - ! Since, in Fortran, there is no null reference, passing the scotchfgraphbuild routine - ! a reference equal to verttab in the velotab or vlbltab fields makes them be considered as missing arrays. - ! The same holds for edlotab when it is passed a reference equal to edgetab. - ! Setting vendtab to refer to one cell after verttab yields the same result, - ! as it is the exact semantics of a compact vertex array. - - call scotchfgraphbuild (scotchgraph(1), 1, nspec_p, & - xadj_tmp(1), xadj_tmp(1), & - elmnts_load_tmp(1), xadj_tmp(1), & - nb_edges_tmp, adjncy_tmp(1), & - adjncy_tmp(1), ier) - if (ier /= 0) stop 'ERROR : MAIN : Cannot build graph' - - ! checks scotch graph - call scotchfgraphcheck (scotchgraph(1), ier) - if (ier /= 0) stop 'Scotch ERROR : MAIN : Invalid check' - - ! computes partition based on scotch graph - if (BALANCE_P_LEVELS_PARTIAL) then - ! partitions level using only a subset of nparts processes - nparts_partial = int( nspec_p / P_LEVEL_PARTIAL_SUBSET_MINIMUM ) - if (nparts_partial > nparts) nparts_partial = nparts - if (nparts_partial < 1) nparts_partial = 1 - - ! user output - print *,' partial partitioning: nparts_partial =',nparts_partial, & - 'using subset maximum',P_LEVEL_PARTIAL_SUBSET_MINIMUM - - ! partitions among subset nparts_partial processes - call scotchfgraphpart (scotchgraph(1), nparts_partial, scotchstrat(1),part_tmp(1),ier) - if (ier /= 0) stop 'Scotch ERROR : MAIN : Cannot part graph' - - else - - ! partitions among all nparts processes - call scotchfgraphpart (scotchgraph(1), nparts, scotchstrat(1),part_tmp(1),ier) - if (ier /= 0) stop 'Scotch ERROR : MAIN : Cannot part graph' + ! computes partition based on scotch graph + if (BALANCE_P_LEVELS_PARTIAL) then + ! partitions level using only a subset of nparts processes + nparts_partial = int( nspec_p / P_LEVEL_PARTIAL_SUBSET_MINIMUM ) + if (nparts_partial > nparts) nparts_partial = nparts + if (nparts_partial < 1) nparts_partial = 1 - endif + ! user output + print *,' partial partitioning: nparts_partial =',nparts_partial, & + 'using subset maximum',P_LEVEL_PARTIAL_SUBSET_MINIMUM - ! frees scotch graph for subsequent calls of scotch again - call scotchfgraphfree (scotchgraph(1),ier) - if (ier /= 0) stop 'Scotch ERROR : MAIN : Cannot free graph' + ! partitions among subset nparts_partial processes + call scotchfgraphpart (scotchgraph(1), nparts_partial, scotchstrat(1),part_tmp(1),ier) + if (ier /= 0) stop 'Scotch ERROR : MAIN : Cannot part graph' - ! stitch partitioning together - do inum = 1,nspec_p - ispec = ispec_global(inum) - part(ispec) = part_tmp(inum) - enddo + else - ! user output: counts number of p-elements in each partition - do iproc = 0,nparts-1 - inum = 0 - do ispec = 1,nspec - if (ispec_p_refine(ispec) == p) then - if (part(ispec) == iproc) inum = inum + 1 - endif - enddo - ! user output - print *,' partition ',iproc, ' has ',inum,' p-elements' - enddo + ! partitions among all nparts processes + call scotchfgraphpart (scotchgraph(1), nparts, scotchstrat(1),part_tmp(1),ier) + if (ier /= 0) stop 'Scotch ERROR : MAIN : Cannot part graph' - enddo ! num_p_level - print * + endif - ! reassemble partitions trying to optimize for communication costs - call remap_partitions(part) + ! frees scotch graph for subsequent calls of scotch again + call scotchfgraphfree (scotchgraph(1),ier) + if (ier /= 0) stop 'Scotch ERROR : MAIN : Cannot free graph' - ! frees memory - deallocate(part_tmp,ispec_global,ispec_local,elmnts_load_tmp,xadj_tmp,adjncy_tmp) + ! stitch partitioning together + do inum = 1,nspec_p + ispec = ispec_global(inum) + part(ispec) = part_tmp(inum) + enddo - end subroutine lts_partition_scotch + ! user output: counts number of p-elements in each partition + do iproc = 0,nparts-1 + inum = 0 + do ispec = 1,nspec + if (ispec_p_refine(ispec) == p) then + if (part(ispec) == iproc) inum = inum + 1 + endif + enddo + ! user output + print *,' partition ',iproc, ' has ',inum,' p-elements' + enddo -#else + enddo ! num_p_level + print * - ! compiled without support - print *,'Scotch Error: Please re-compile with -DUSE_SCOTCH to use SCOTCH partitioning.' - stop 'No SCOTCH support' + ! reassemble partitions trying to optimize for communication costs + call remap_partitions(part) -#endif + ! frees memory + deallocate(part_tmp,ispec_global,ispec_local,elmnts_load_tmp,xadj_tmp,adjncy_tmp) - end subroutine partition_scotch + end subroutine lts_partition_scotch +#endif ! !---------------------------------------------------------------------------------------------- diff --git a/src/decompose_mesh/read_mesh_files.F90 b/src/decompose_mesh/read_mesh_files.F90 index 5051991ad..a15a2abc2 100644 --- a/src/decompose_mesh/read_mesh_files.F90 +++ b/src/decompose_mesh/read_mesh_files.F90 @@ -277,7 +277,7 @@ subroutine read_mesh_files() if (len_trim(line) == 0) cycle if (line(1:1) == '#' .or. line(1:1) == '!') cycle - read(line,*,iostat=ier) idomain_id,num_mat + read(line,*) idomain_id,num_mat enddo if (ier /= 0) stop 'Error reading in defined materials in nummaterial_velocity_file' @@ -322,7 +322,7 @@ subroutine read_mesh_files() ! reads poroelastic file line do while (ier == 0) - read(98,'(A)',iostat=ier) line + read(97,'(A)',iostat=ier) line if (ier /= 0) exit ! skip empty/comment lines @@ -333,9 +333,10 @@ subroutine read_mesh_files() read(line,*,iostat=ier) rhos,rhof,phi,tort,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,eta,mufr if (ier /= 0) then stop 'Error reading nummaterial_poroelastic_file, please check if it has the required format...' - else - exit endif + + ! only read 1 valid line at a time + if (ier == 0) exit enddo if (ier /= 0) stop 'Error reading in materials in nummaterial_poroelastic_file' diff --git a/src/decompose_mesh/rules.mk b/src/decompose_mesh/rules.mk index 33671197e..1f43585cf 100644 --- a/src/decompose_mesh/rules.mk +++ b/src/decompose_mesh/rules.mk @@ -53,6 +53,7 @@ decompose_mesh_MODULES = \ $(FC_MODDIR)/module_database.$(FC_MODEXT) \ $(FC_MODDIR)/module_mesh.$(FC_MODEXT) \ $(FC_MODDIR)/module_partition.$(FC_MODEXT) \ + $(FC_MODDIR)/scotch_par.$(FC_MODEXT) \ $(EMPTY_MACRO) decompose_mesh_SHARED_OBJECTS = \ diff --git a/src/generate_databases/calc_jacobian.f90 b/src/generate_databases/calc_jacobian.f90 index 0c1c5aff7..c4bd1127d 100644 --- a/src/generate_databases/calc_jacobian.f90 +++ b/src/generate_databases/calc_jacobian.f90 @@ -31,7 +31,7 @@ subroutine calc_jacobian(myrank,xix_elem,xiy_elem,xiz_elem, & gammax_elem,gammay_elem,gammaz_elem,jacobian_elem, & xelm,yelm,zelm,dershape3D) - use generate_databases_par, only: NGNOD,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,SIZE_REAL,ZERO + use generate_databases_par, only: NGNOD,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,ZERO implicit none @@ -168,7 +168,7 @@ end subroutine calc_coords subroutine check_element_regularity(xelm,yelm,zelm,any_regular_elem,cube_edge_size_squared, & nspec_irregular,ispec,nspec,irregular_element_number,ANY_FAULT_IN_THIS_PROC) - use generate_databases_par, only: NGNOD,CUSTOM_REAL,USE_MESH_COLORING_GPU + use generate_databases_par, only: NGNOD,USE_MESH_COLORING_GPU real, dimension(NGNOD) :: xelm,yelm,zelm logical :: any_regular_elem,ANY_FAULT_IN_THIS_PROC @@ -178,12 +178,18 @@ subroutine check_element_regularity(xelm,yelm,zelm,any_regular_elem,cube_edge_si !local real :: dist1_sq,dist2_sq,dist3_sq + double precision :: threshold + double precision,parameter :: threshold_percentage = 1.e-5 !checks if the potential cube has the same size as the previous ones dist1_sq = (xelm(2)-xelm(1))**2 + (yelm(2)-yelm(1))**2 +(zelm(2)-zelm(1))**2 - if (NGNOD == 27 .or. ANY_FAULT_IN_THIS_PROC .or. USE_MESH_COLORING_GPU .or. & - (any_regular_elem .and. ( abs(dist1_sq - cube_edge_size_squared) > (1e-5)*cube_edge_size_squared ))) then + ! sets irregular element (default for NGNOD 27 and others) + threshold = threshold_percentage*cube_edge_size_squared + if (NGNOD == 27 & + .or. ANY_FAULT_IN_THIS_PROC & + .or. USE_MESH_COLORING_GPU & + .or. (any_regular_elem .and. ( abs(dist1_sq - cube_edge_size_squared) > threshold ))) then irregular_element_number(ispec) = ispec - (nspec - nspec_irregular) return endif @@ -199,8 +205,9 @@ subroutine check_element_regularity(xelm,yelm,zelm,any_regular_elem,cube_edge_si dist2_sq = (xelm(5)-xelm(1))**2 + (yelm(5)-yelm(1))**2 +(zelm(5)-zelm(1))**2 dist3_sq = (xelm(4)-xelm(1))**2 + (yelm(4)-yelm(1))**2 +(zelm(4)-zelm(1))**2 - if (abs(dist2_sq - dist1_sq) < 1e-5*dist1_sq .and. abs(dist3_sq - dist1_sq) < 1e-5*dist1_sq) then + threshold = threshold_percentage*dist1_sq + if (abs(dist2_sq - dist1_sq) < threshold .and. abs(dist3_sq - dist1_sq) < threshold) then ! test if first cube found in mesh if (.not. any_regular_elem ) then cube_edge_size_squared = dist1_sq diff --git a/src/generate_databases/create_mass_matrices.f90 b/src/generate_databases/create_mass_matrices.f90 index 1fa912838..8ae8bfa83 100644 --- a/src/generate_databases/create_mass_matrices.f90 +++ b/src/generate_databases/create_mass_matrices.f90 @@ -29,7 +29,7 @@ subroutine create_mass_matrices(nglob,nspec,ibool,PML_CONDITIONS,STACEY_ABSORBIN ! returns precomputed mass matrix in rmass array - use constants, only: NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,SIZE_REAL + use constants, only: NGLLX,NGLLY,NGLLZ,CUSTOM_REAL use create_regions_mesh_ext_par @@ -185,7 +185,7 @@ subroutine create_mass_matrices_ocean_load(nglob,nspec,ibool) use generate_databases_par, only: & APPROXIMATE_OCEAN_LOAD,TOPOGRAPHY, & NX_TOPO,NY_TOPO,itopo_bathy,myrank, & - NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,SIZE_REAL,NGLLSQUARE,IMAIN, & + NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,NGLLSQUARE,IMAIN, & MINIMUM_THICKNESS_3D_OCEANS,RHO_APPROXIMATE_OCEAN_LOAD use create_regions_mesh_ext_par @@ -303,7 +303,7 @@ subroutine create_mass_matrices_Stacey(nglob,nspec,ibool) ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix - use generate_databases_par, only: DT,NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,SIZE_REAL,NGLLSQUARE + use generate_databases_par, only: DT,NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,NGLLSQUARE use create_regions_mesh_ext_par implicit none @@ -421,7 +421,7 @@ subroutine create_mass_matrices_pml_elastic(nspec,ibool) use generate_databases_par, only: DT,is_CPML,CPML_regions,d_store_x,d_store_y,d_store_z, & K_store_x,K_store_y,K_store_z,nspec_cpml,CPML_to_spec, & - NGLLX,NGLLY,NGLLZ,SIZE_REAL,CUSTOM_REAL, & + NGLLX,NGLLY,NGLLZ,CUSTOM_REAL, & CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, & CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ @@ -640,7 +640,7 @@ subroutine create_mass_matrices_pml_acoustic(nspec,ibool) use generate_databases_par, only: DT,is_CPML,CPML_regions,d_store_x,d_store_y,d_store_z, & K_store_x,K_store_y,K_store_z,nspec_cpml,CPML_to_spec, & - NGLLX,NGLLY,NGLLZ,SIZE_REAL,CUSTOM_REAL, & + NGLLX,NGLLY,NGLLZ,CUSTOM_REAL, & CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, & CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ diff --git a/src/generate_databases/create_regions_mesh.f90 b/src/generate_databases/create_regions_mesh.f90 index 549c0ad97..8b9da2f5f 100644 --- a/src/generate_databases/create_regions_mesh.f90 +++ b/src/generate_databases/create_regions_mesh.f90 @@ -1021,7 +1021,7 @@ subroutine crm_setup_moho( myrank,nspec, & nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & nodes_coords_ext_mesh,nnodes_ext_mesh,ibool ) - use generate_databases_par, only: NGNOD2D,NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,SIZE_REAL,IMAIN, & + use generate_databases_par, only: NGNOD2D,NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,IMAIN, & NDIM,NGLLSQUARE,NGNOD2D_FOUR_CORNERS use create_regions_mesh_ext_par diff --git a/src/generate_databases/fault_generate_databases.f90 b/src/generate_databases/fault_generate_databases.f90 index f0463c821..e42aaeaef 100644 --- a/src/generate_databases/fault_generate_databases.f90 +++ b/src/generate_databases/fault_generate_databases.f90 @@ -105,9 +105,17 @@ subroutine fault_read_input(prname,myrank) open(unit=IIN_PAR,file=IN_DATA_FILES(1:len_trim(IN_DATA_FILES))//'Par_file_faults',status='old',action='read',iostat=ier) if (ier == 0) then read(IIN_PAR,*) nb - if (myrank == 0) write(IMAIN,*) ' ... reading ', nb,' faults from file DATA/Par_file_faults' + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) ' ... reading ', nb,' faults from file DATA/Par_file_faults' + write(IMAIN,*) + endif else - if (myrank == 0) write(IMAIN,*) 'File DATA/Par_file_faults not found: assuming that there are no faults' + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) 'File DATA/Par_file_faults not found: assuming that there are no faults' + write(IMAIN,*) + endif close(IIN_PAR) endif diff --git a/src/generate_databases/get_model.F90 b/src/generate_databases/get_model.F90 index e53988966..a8245171f 100644 --- a/src/generate_databases/get_model.F90 +++ b/src/generate_databases/get_model.F90 @@ -28,8 +28,7 @@ subroutine get_model(myrank) use generate_databases_par, only: IMODEL, & - IMODEL_DEFAULT,IMODEL_GLL,IMODEL_1D_PREM,IMODEL_1D_CASCADIA,IMODEL_1D_SOCAL, & - IMODEL_SALTON_TROUGH,IMODEL_TOMO,IMODEL_USER_EXTERNAL,IMODEL_IPATI,IMODEL_IPATI_WATER, & + IMODEL_SALTON_TROUGH,IMODEL_TOMO,IMODEL_USER_EXTERNAL, & IMODEL_COUPLED, & IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC, & nspec => NSPEC_AB,ibool,mat_ext_mesh, & diff --git a/src/generate_databases/model_coupled.f90 b/src/generate_databases/model_coupled.f90 index 7a87dc774..6417a7967 100644 --- a/src/generate_databases/model_coupled.f90 +++ b/src/generate_databases/model_coupled.f90 @@ -232,7 +232,7 @@ subroutine model_1D_coupling(x_eval,y_eval,z_eval,rho_final,vp_final,vs_final,r1 double precision :: r1,radius double precision :: rho,vp,vs - double precision, parameter :: Xtol = 1d-2 + !double precision, parameter :: Xtol = 1d-2 radius = dsqrt(x_eval**2 + y_eval**2 + (z_eval+ZREF)**2) radius = radius / 1000.d0 diff --git a/src/generate_databases/model_sep.f90 b/src/generate_databases/model_sep.f90 index 21a6391b2..2f642ce86 100644 --- a/src/generate_databases/model_sep.f90 +++ b/src/generate_databases/model_sep.f90 @@ -40,11 +40,13 @@ module model_sep_mod subroutine model_sep() - use generate_databases_par, only: NGLLX, NGLLY, NGLLZ, & - SEP_MODEL_DIRECTORY, FOUR_THIRDS, myrank, IMAIN - use create_regions_mesh_ext_par, only: rhostore, rho_vp, rho_vs, & - kappastore, mustore, & - rho_vpI, rho_vsI, rhoarraystore + use generate_databases_par, only: & + SEP_MODEL_DIRECTORY, FOUR_THIRDS, myrank, IMAIN + + use create_regions_mesh_ext_par, only: & + rhostore, rho_vp, rho_vs, & + kappastore, mustore, & + rho_vpI, rho_vsI, rhoarraystore implicit none real(kind=4), allocatable, dimension(:,:,:) :: vp_sep, vs_sep, rho_sep diff --git a/src/generate_databases/pml_set_local_dampingcoeff.f90 b/src/generate_databases/pml_set_local_dampingcoeff.f90 index 1874be338..94f6d6155 100644 --- a/src/generate_databases/pml_set_local_dampingcoeff.f90 +++ b/src/generate_databases/pml_set_local_dampingcoeff.f90 @@ -33,8 +33,8 @@ subroutine pml_set_local_dampingcoeff(myrank,xstore,ystore,zstore) use generate_databases_par, only: ibool,NGLOB_AB,d_store_x,d_store_y,d_store_z, & K_store_x,K_store_y,K_store_z,alpha_store_x,alpha_store_y,alpha_store_z,CPML_to_spec, & - CPML_width_x,CPML_width_y,CPML_width_z,min_distance_between_CPML_parameter,NPOWER, & - CUSTOM_REAL,SIZE_REAL,NGLLX,NGLLY,NGLLZ,nspec_cpml,PML_INSTEAD_OF_FREE_SURFACE, & + CPML_width_x,CPML_width_y,CPML_width_z,min_distance_between_CPML_parameter, & + CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,nspec_cpml,PML_INSTEAD_OF_FREE_SURFACE, & IMAIN,CPML_REGIONS,f0_FOR_PML,PI, & CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ, & SIMULATION_TYPE,SAVE_FORWARD,nspec => NSPEC_AB,is_CPML, & diff --git a/src/generate_databases/save_arrays_solver.F90 b/src/generate_databases/save_arrays_solver.F90 index f76a9e4c0..91c528baf 100644 --- a/src/generate_databases/save_arrays_solver.F90 +++ b/src/generate_databases/save_arrays_solver.F90 @@ -32,7 +32,7 @@ subroutine save_arrays_solver_ext_mesh(nspec,nglob,APPROXIMATE_OCEAN_LOAD,ibool, max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, & SAVE_MESH_FILES,ANISOTROPY) - use generate_databases_par, only: NGLLX,NGLLY,NGLLZ,NGLLSQUARE,IMAIN,IOUT, & + use generate_databases_par, only: NGLLX,NGLLY,NGLLZ,IOUT, & nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, & ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & SIMULATION_TYPE,SAVE_FORWARD,mask_ibool_interior_domain, & @@ -726,7 +726,7 @@ end subroutine save_arrays_solver_files subroutine save_arrays_solver_injection_boundary(nspec,ibool) - use generate_databases_par, only: myrank,NGLLX,NGLLY,NGLLZ,NGLLSQUARE,IMAIN,IOUT,FOUR_THIRDS + use generate_databases_par, only: myrank,NGLLX,NGLLY,NGLLZ,NGLLSQUARE,IMAIN,IOUT use create_regions_mesh_ext_par diff --git a/src/generate_databases/setup_color_perm.f90 b/src/generate_databases/setup_color_perm.f90 index 705e1af83..a5e47f75c 100644 --- a/src/generate_databases/setup_color_perm.f90 +++ b/src/generate_databases/setup_color_perm.f90 @@ -132,7 +132,7 @@ subroutine setup_color(myrank,nspec,nglob,ibool,perm,ispec_is_d,idomain, & ! sets up mesh coloring - use generate_databases_par, only: NGLLX,NGLLY,NGLLZ,IMAIN,USE_MESH_COLORING_GPU,MAX_NUMBER_OF_COLORS + use generate_databases_par, only: NGLLX,NGLLY,NGLLZ,IMAIN,MAX_NUMBER_OF_COLORS use create_regions_mesh_ext_par implicit none diff --git a/src/inverse_problem_for_model/input_output/IO_model_mod.f90 b/src/inverse_problem_for_model/input_output/IO_model_mod.f90 index a512243bf..ed8d0c6d5 100644 --- a/src/inverse_problem_for_model/input_output/IO_model_mod.f90 +++ b/src/inverse_problem_for_model/input_output/IO_model_mod.f90 @@ -611,7 +611,7 @@ end subroutine DumpModelArray subroutine import_FD_model_ACOUSTIC(fd_grid) use specfem_par, only: myrank, xstore, ystore, zstore, kappastore, rhostore, ibool, & - NGLLX, NGLLY, NGLLZ, NSPEC_AB, FOUR_THIRDS, MAX_STRING_LEN + NGLLX, NGLLY, NGLLZ, NSPEC_AB, MAX_STRING_LEN implicit none @@ -875,7 +875,7 @@ subroutine import_FD_model_ANISO(fd_grid, inversion_param) use constants, only: MAX_STRING_LEN use specfem_par, only: myrank, xstore, ystore, zstore, rhostore, ibool, & - NGLLX, NGLLY, NGLLZ, NSPEC_AB, FOUR_THIRDS, MAX_STRING_LEN + NGLLX, NGLLY, NGLLZ, NSPEC_AB, MAX_STRING_LEN !use specfem_par_elastic, only: rho_vp, rho_vs here i commented because they are already called in head of module use interpolation_mod, only: trilin_interp diff --git a/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 b/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 index 0c1922683..f87e31675 100644 --- a/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 +++ b/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 @@ -495,10 +495,11 @@ end subroutine compute_source_coeff ! compute mass matrix !--------------------------------------------------------------- subroutine create_mass_matrices_Stacey_duplication_routine() + use constants, only: NGLLSQUARE use shared_parameters, only: NPROC, DT, STACEY_ABSORBING_CONDITIONS, PML_CONDITIONS use specfem_par, only: kappastore, rhostore, ibool, wxgll, wygll, wzgll, jacobian, & - NGLLX, NGLLY, NGLLZ, NSPEC_AB, NGLOB_AB, FOUR_THIRDS, num_abs_boundary_faces, abs_boundary_ispec, abs_boundary_ijk, & + NGLLX, NGLLY, NGLLZ, NSPEC_AB, NGLOB_AB, num_abs_boundary_faces, abs_boundary_ispec, abs_boundary_ijk, & abs_boundary_normal, abs_boundary_jacobian2Dw, num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, my_neighbors_ext_mesh use specfem_par_elastic, only: rmass, rmassx, rmassy, rmassz, rho_vp, rho_vs, ispec_is_elastic, ELASTIC_SIMULATION diff --git a/src/meshfem3D/check_mesh_quality.f90 b/src/meshfem3D/check_mesh_quality.f90 index 190c353ab..118fcdbc6 100644 --- a/src/meshfem3D/check_mesh_quality.f90 +++ b/src/meshfem3D/check_mesh_quality.f90 @@ -405,11 +405,10 @@ subroutine create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NGLOB,VP_MAX,dt_s stability,distmin,distmax) use constants + use constants_meshfem3D, only: NGLLX_M implicit none - include "constants_meshfem3D.h" - integer :: NSPEC,NGLOB double precision, dimension(NGLOB) :: x,y,z diff --git a/src/meshfem3D/compute_parameters.f90 b/src/meshfem3D/compute_parameters.f90 index b424b66a1..d2b431048 100644 --- a/src/meshfem3D/compute_parameters.f90 +++ b/src/meshfem3D/compute_parameters.f90 @@ -34,11 +34,10 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, & USE_REGULAR_MESH,NDOUBLINGS,ner_doublings) use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M implicit none - include "constants_meshfem3D.h" - ! parameters read from parameter file integer,intent(in) :: NER integer,intent(in) :: NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA diff --git a/src/meshfem3D/create_interfaces_mesh.f90 b/src/meshfem3D/create_interfaces_mesh.f90 index 5897687cb..0c88278d4 100644 --- a/src/meshfem3D/create_interfaces_mesh.f90 +++ b/src/meshfem3D/create_interfaces_mesh.f90 @@ -460,7 +460,7 @@ subroutine get_interfaces_mesh_count() number_of_interfaces,max_npx_interface,max_npy_interface, & number_of_layers,ner_layer - use constants, only: IMAIN,IIN,MF_IN_DATA_FILES,DONT_IGNORE_JUNK,IUTM2LONGLAT,MAX_STRING_LEN + use constants, only: IMAIN,IIN,MF_IN_DATA_FILES,DONT_IGNORE_JUNK,MAX_STRING_LEN implicit none diff --git a/src/meshfem3D/create_meshfem_mesh.f90 b/src/meshfem3D/create_meshfem_mesh.f90 index 01c89ffe7..ceada7e9a 100644 --- a/src/meshfem3D/create_meshfem_mesh.f90 +++ b/src/meshfem3D/create_meshfem_mesh.f90 @@ -43,6 +43,14 @@ module create_meshfem_par ! MPI cut-planes parameters along xi and along eta logical, dimension(:,:), allocatable :: iMPIcut_xi,iMPIcut_eta + ! flag indicating whether point is in the sediments + ! logical point_is_in_sediments + logical, dimension(:,:,:,:), allocatable :: flag_sediments + logical, dimension(:), allocatable :: not_fully_in_bedrock + + ! boundary parameters locator + integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top + end module create_meshfem_par ! @@ -52,76 +60,38 @@ end module create_meshfem_par subroutine create_meshfem_mesh() - use create_meshfem_par, only: nodes_coords,ispec_material_id,iboun,iMPIcut_xi,iMPIcut_eta +! create the different regions of the mesh - use meshfem3D_par, only: NSPEC_AB,xgrid,ygrid,zgrid,ibool, & + use constants, only: IMAIN + + use meshfem3D_par, only: & + ibool, & xstore,ystore,zstore, & iproc_xi_current,iproc_eta_current,addressing,nspec, & - NGLOB_AB,npointot, & - NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NER, & NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, & NPROC_XI,NPROC_ETA, & - nsubregions,subregions,NMATERIALS,material_properties, & + NMATERIALS,material_properties, & myrank, sizeprocs, prname, & - LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,NEX_XI,NEX_ETA, & + LOCAL_PATH, & CREATE_ABAQUS_FILES,CREATE_DX_FILES,CREATE_VTK_FILES, & - USE_REGULAR_MESH,NDOUBLINGS,ner_doublings, & ADIOS_ENABLED, ADIOS_FOR_DATABASES, & nspec_CPML,is_CPML,CPML_to_spec,CPML_regions - ! create the different regions of the mesh - use constants, only: MF_IN_DATA_FILES,MAX_STRING_LEN,NGNOD_EIGHT_CORNERS, & - IMAIN,IOVTK,CUSTOM_REAL,HUGEVAL,TINYVAL,NDIM - - use constants_meshfem3D, only: NSPEC_DOUBLING_SUPERBRICK,NGLOB_DOUBLING_SUPERBRICK, & - IFLAG_BASEMENT_TOPO,IFLAG_ONE_LAYER_TOPOGRAPHY, & - NGLLCUBE_M,NGLLX_M,NGLLY_M,NGLLZ_M + use create_meshfem_par use adios_manager_mod, only: adios_setup,adios_cleanup implicit none ! local parameters - ! auxiliary variables to generate the mesh - integer :: isubregion,doubling_index,nmeshregions - integer, dimension(:,:,:), allocatable :: material_num - - ! topology of the elements - integer :: iaddx(NGNOD_EIGHT_CORNERS) - integer :: iaddy(NGNOD_EIGHT_CORNERS) - integer :: iaddz(NGNOD_EIGHT_CORNERS) - - double precision :: xelm(NGNOD_EIGHT_CORNERS) - double precision :: yelm(NGNOD_EIGHT_CORNERS) - double precision :: zelm(NGNOD_EIGHT_CORNERS) - - ! variables for creating array ibool (some arrays also used for AVS or DX files) - integer, dimension(:), allocatable :: iglob,locval - logical, dimension(:), allocatable :: ifseg - double precision, dimension(:), allocatable :: xp,yp,zp - integer :: nglob - ! boundary parameters locator - integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top - ! number of elements on the boundaries integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax - ! flag indicating whether point is in the sediments - ! logical point_is_in_sediments - logical, dimension(:,:,:,:), allocatable :: flag_sediments - logical, dimension(:), allocatable :: not_fully_in_bedrock - ! material properties double precision :: VP_MAX - ! doublings zone - integer :: nspec_sb - logical, dimension(NSPEC_DOUBLING_SUPERBRICK,6) :: iboun_sb - integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick - double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick - ! ************** !--- Initialize ADIOS and setup the buffer size @@ -139,7 +109,7 @@ subroutine create_meshfem_mesh() call cmm_create_mesh_elements() ! creates indirect addressing - call cmm_create_addressing() + call cmm_create_addressing(nglob) ! determine cavity call cmm_determine_cavity(nglob) @@ -205,248 +175,311 @@ subroutine create_meshfem_mesh() ! frees memory deallocate(ispec_material_id) - deallocate(locval,ifseg,xp,yp,zp) deallocate(nodes_coords) deallocate(ibool,xstore,ystore,zstore) deallocate(is_CPML) deallocate(CPML_to_spec) deallocate(CPML_regions) -contains + end subroutine create_meshfem_mesh ! !---------------------------------------------------------------------------------------------------- ! - subroutine cmm_allocate_arrays() + subroutine cmm_allocate_arrays() - implicit none - ! local parameters - integer :: ier + use constants, only: IMAIN + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NGLLCUBE_M - ! user output - if (myrank == 0) then - write(IMAIN,*) 'allocating mesh arrays' - write(IMAIN,*) - call flush_IMAIN() - endif + use meshfem3D_par, only: NSPEC_AB,nspec,ibool, & + xstore,ystore,zstore,npointot,myrank, & + NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP + + use create_meshfem_par + + implicit none + ! local parameters + integer :: ier + + ! user output + if (myrank == 0) then + write(IMAIN,*) 'allocating mesh arrays' + write(IMAIN,*) + call flush_IMAIN() + endif + + ! assign theoretical number of elements + nspec = NSPEC_AB + + ! compute maximum number of points + npointot = nspec * NGLLCUBE_M + + ! make sure everybody is synchronized + call synchronize_all() - ! assign theoretical number of elements - nspec = NSPEC_AB - - ! compute maximum number of points - npointot = nspec * NGLLCUBE_M - - ! make sure everybody is synchronized - call synchronize_all() - - ! use dynamic allocation to allocate memory for arrays - allocate(ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1257') - if (ier /= 0) stop 'Error allocating array ibool' - - allocate(xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1258') - if (ier /= 0) stop 'Error allocating array xstore' - allocate(ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1259') - if (ier /= 0) stop 'Error allocating array ystore' - allocate(zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1260') - ! exit if there is not enough memory to allocate all the arrays - if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - - ! flag indicating whether point is in the sediments - allocate(flag_sediments(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1261') - if (ier /= 0) stop 'Error allocating array flag_sediments' - allocate(not_fully_in_bedrock(nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1262') - if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - - ! boundary locator - allocate(iboun(6,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1263') - if (ier /= 0) stop 'Error allocating array iboun' - - ! boundary parameters locator - allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1264') - if (ier /= 0) stop 'Error allocating array ibelm_xmin' - allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1265') - if (ier /= 0) stop 'Error allocating array ibelm_xmax' - allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1266') - if (ier /= 0) stop 'Error allocating array ibelm_ymin' - allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1267') - if (ier /= 0) stop 'Error allocating array ibelm_ymax' - allocate(ibelm_bottom(NSPEC2D_BOTTOM),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1268') - if (ier /= 0) stop 'Error allocating array ibelm_bottom' - allocate(ibelm_top(NSPEC2D_TOP),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1269') - if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - - ! MPI cut-planes parameters along xi and along eta - allocate(iMPIcut_xi(2,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1270') - if (ier /= 0) stop 'Error allocating array iMPIcut_xi' - allocate(iMPIcut_eta(2,nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1271') - if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - - ! allocate memory for arrays - allocate(iglob(npointot),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1272') - if (ier /= 0) stop 'Error allocating array iglob' - allocate(locval(npointot),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1273') - if (ier /= 0) stop 'Error allocating array locval' - allocate(ifseg(npointot),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1274') - if (ier /= 0) stop 'Error allocating array ifseg' - allocate(xp(npointot),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1275') - if (ier /= 0) stop 'Error allocating array xp' - allocate(yp(npointot),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1276') - if (ier /= 0) stop 'Error allocating array yp' - allocate(zp(npointot),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1277') - if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - - ! allocate material ids array - allocate(material_num(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1278') - if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - - allocate(ispec_material_id(nspec),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1279') - if (ier /= 0) stop 'Error allocating array ispec_material_id' - - ! synchronize - call synchronize_all() - - end subroutine cmm_allocate_arrays + ! use dynamic allocation to allocate memory for arrays + allocate(ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1257') + + allocate(xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1258') + allocate(ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1259') + allocate(zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1260') + + ! flag indicating whether point is in the sediments + allocate(flag_sediments(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1261') + allocate(not_fully_in_bedrock(nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1262') + + ! boundary locator + allocate(iboun(6,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1263') + + ! boundary parameters locator + allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1264') + allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1265') + allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1266') + allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1267') + allocate(ibelm_bottom(NSPEC2D_BOTTOM),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1268') + allocate(ibelm_top(NSPEC2D_TOP),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1269') + + ! MPI cut-planes parameters along xi and along eta + allocate(iMPIcut_xi(2,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1270') + allocate(iMPIcut_eta(2,nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1271') + + allocate(ispec_material_id(nspec),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1279') + + ! synchronize + call synchronize_all() + + end subroutine cmm_allocate_arrays ! !---------------------------------------------------------------------------------------------------- ! - subroutine cmm_create_mesh_elements() - - implicit none - ! local parameters - integer :: ix,iy,ir,ir1,ir2,dir - integer :: ix1,ix2,dix,iy1,iy2,diy - integer :: iax,iay,iar,ia - integer :: imaterial_number - integer :: ioffset_x,ioffset_y,ioffset_z - integer :: ispec,ispec_superbrick - double precision :: x_min,x_max,y_min,y_max,z_min,z_max - double precision :: x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob - real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max,elemsize_min_glob,elemsize_max_glob - real(kind=CUSTOM_REAL) :: dist,x1,y1,z1,x2,y2,z2 - integer :: i1,i2,j1,j2,k1,k2,i,j,k - - ! generate the elements in all the regions of the mesh - ispec = 0 - - ! to check that we assign a material to each element - material_num(:,:,:) = -1000 ! dummy value - ispec_material_id(:) = -1000 + subroutine cmm_create_mesh_elements() + + use constants, only: CUSTOM_REAL,IMAIN,NGNOD_EIGHT_CORNERS,HUGEVAL + + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M, & + NGLOB_DOUBLING_SUPERBRICK,NSPEC_DOUBLING_SUPERBRICK, & + IFLAG_BASEMENT_TOPO,IFLAG_ONE_LAYER_TOPOGRAPHY + + use meshfem3D_par, only: myrank,xstore,ystore,zstore, & + xgrid,ygrid,zgrid, & + UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, & + NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, & + NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NER, & + iproc_xi_current,iproc_eta_current, & + nspec, & + USE_REGULAR_MESH,NDOUBLINGS,ner_doublings, & + nsubregions,subregions + + use create_meshfem_par + + implicit none + ! local parameters + integer :: ix,iy,ir,ir1,ir2,dir + integer :: ix1,ix2,dix,iy1,iy2,diy + integer :: iax,iay,iar,ia + integer :: imaterial_number + integer :: ioffset_x,ioffset_y,ioffset_z + integer :: ispec,ispec_superbrick + integer :: i1,i2,j1,j2,k1,k2,i,j,k,ier + + double precision :: x_min,x_max,y_min,y_max,z_min,z_max + double precision :: x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob + + real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max,elemsize_min_glob,elemsize_max_glob + real(kind=CUSTOM_REAL) :: dist,x1,y1,z1,x2,y2,z2 + + ! auxiliary variables to generate the mesh + integer :: isubregion,doubling_index,nmeshregions + + ! topology of the elements + integer :: iaddx(NGNOD_EIGHT_CORNERS) + integer :: iaddy(NGNOD_EIGHT_CORNERS) + integer :: iaddz(NGNOD_EIGHT_CORNERS) + double precision :: xelm(NGNOD_EIGHT_CORNERS) + double precision :: yelm(NGNOD_EIGHT_CORNERS) + double precision :: zelm(NGNOD_EIGHT_CORNERS) + + integer, dimension(:,:,:), allocatable :: material_num + + ! doublings zone + integer :: nspec_sb + logical, dimension(NSPEC_DOUBLING_SUPERBRICK,6) :: iboun_sb + integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick + double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick + + ! allocate material ids array + allocate(material_num(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1278') + + ! generate the elements in all the regions of the mesh + ispec = 0 + + ! to check that we assign a material to each element + material_num(:,:,:) = -1000 ! dummy value + ispec_material_id(:) = -1000 + + ! user output + if (myrank == 0) then + write(IMAIN,*) 'number of subregions = ',nsubregions + call flush_IMAIN() + endif + + do isubregion = 1,nsubregions ! user output if (myrank == 0) then - write(IMAIN,*) 'number of subregions = ',nsubregions + write(IMAIN,*) ' defining subregion ',isubregion call flush_IMAIN() endif - do isubregion = 1,nsubregions - ! user output - if (myrank == 0) then - write(IMAIN,*) ' defining subregion ',isubregion - call flush_IMAIN() - endif - - call define_model_regions(NEX_PER_PROC_XI,NEX_PER_PROC_ETA,iproc_xi_current,iproc_eta_current, & - isubregion,nsubregions,subregions, & - iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, & - imaterial_number) + call define_model_regions(NEX_PER_PROC_XI,NEX_PER_PROC_ETA,iproc_xi_current,iproc_eta_current, & + isubregion,nsubregions,subregions, & + iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, & + imaterial_number) - ! loop on all the mesh points in current subregion - do ir = ir1,ir2,dir - do iy = iy1,iy2,diy - do ix = ix1,ix2,dix - material_num(ir,ix,iy) = imaterial_number - enddo + ! loop on all the mesh points in current subregion + do ir = ir1,ir2,dir + do iy = iy1,iy2,diy + do ix = ix1,ix2,dix + material_num(ir,ix,iy) = imaterial_number enddo enddo - enddo - if (USE_REGULAR_MESH) then - nmeshregions = 1 - else - call define_superbrick(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb) - nspec_sb = NSPEC_DOUBLING_SUPERBRICK - - ! mesh subregion: - ! example NDOUBLINGS = 1: isubregion = 1,2,3 - ! example NDOUBLINGS = 2: isubregion = 1,2,3,4,5 - ! the values goes from 1 to 2 * NDOUBLINGS + 1 - nmeshregions = 2 * NDOUBLINGS + 1 - - ! hard-coded way - !if (NDOUBLINGS == 1) then - ! nmeshregions = 3 - !else if (NDOUBLINGS == 2) then - ! nmeshregions = 5 - !else - ! stop 'Wrong number of mesh doubling regions' - !endif - endif + enddo + + if (USE_REGULAR_MESH) then + nmeshregions = 1 + else + call define_superbrick(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb) + nspec_sb = NSPEC_DOUBLING_SUPERBRICK + + ! mesh subregion: + ! example NDOUBLINGS = 1: isubregion = 1,2,3 + ! example NDOUBLINGS = 2: isubregion = 1,2,3,4,5 + ! the values goes from 1 to 2 * NDOUBLINGS + 1 + nmeshregions = 2 * NDOUBLINGS + 1 + + ! hard-coded way + !if (NDOUBLINGS == 1) then + ! nmeshregions = 3 + !else if (NDOUBLINGS == 2) then + ! nmeshregions = 5 + !else + ! stop 'Wrong number of mesh doubling regions' + !endif + endif + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) 'number of mesh regions = ',nmeshregions + call flush_IMAIN() + endif + + do isubregion = 1,nmeshregions ! user output if (myrank == 0) then - write(IMAIN,*) - write(IMAIN,*) 'number of mesh regions = ',nmeshregions + if (modulo(isubregion,2) == 1) then + write(IMAIN,*) ' creating mesh region ',isubregion + else + write(IMAIN,*) ' creating mesh region ',isubregion,' with doubling layer' + endif call flush_IMAIN() endif - do isubregion = 1,nmeshregions - ! user output - if (myrank == 0) then - if (modulo(isubregion,2) == 1) then - write(IMAIN,*) ' creating mesh region ',isubregion - else - write(IMAIN,*) ' creating mesh region ',isubregion,' with doubling layer' - endif - call flush_IMAIN() - endif + ! define shape of elements + call define_mesh_regions(myrank,USE_REGULAR_MESH,isubregion,NER,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, & + iproc_xi_current,iproc_eta_current, & + NDOUBLINGS,ner_doublings, & + iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar) - ! define shape of elements - call define_mesh_regions(myrank,USE_REGULAR_MESH,isubregion,NER,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, & - iproc_xi_current,iproc_eta_current, & - NDOUBLINGS,ner_doublings, & - iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar) + ! loop on all the mesh points in current subregion + if (modulo(isubregion,2) == 1) then - ! loop on all the mesh points in current subregion - if (modulo(isubregion,2) == 1) then + ! Regular subregion case + + do ir = ir1,ir2,dir + do iy = iy1,iy2,diy + do ix = ix1,ix2,dix + ! loop over the NGNOD_EIGHT_CORNERS nodes + do ia = 1,NGNOD_EIGHT_CORNERS + ! define topological coordinates of this mesh point + ioffset_x = ix+iax*iaddx(ia) + ioffset_y = iy+iay*iaddy(ia) + ioffset_z = ir+iar*iaddz(ia) + + xelm(ia) = xgrid(ioffset_z,ioffset_x,ioffset_y) + yelm(ia) = ygrid(ioffset_z,ioffset_x,ioffset_y) + zelm(ia) = zgrid(ioffset_z,ioffset_x,ioffset_y) + enddo + + ! add one spectral element to the list and store its material number + ispec = ispec + 1 + + if (ispec > nspec) then + call exit_MPI(myrank,'ispec greater than nspec in mesh creation') + endif + + ! check if element is on topography + if ((ir == ir2) .and. (isubregion == nmeshregions)) then + ! sets index to topography layer + doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY + else + ! sets index to deeper, basement layers + doubling_index = IFLAG_BASEMENT_TOPO + endif + + ispec_material_id(ispec) = material_num(ir,ix,iy) + + ! store coordinates + call store_coords(xstore,ystore,zstore,xelm,yelm,zelm,ispec,nspec) + + ! detect mesh boundaries + call get_flags_boundaries(nspec,iproc_xi_current,iproc_eta_current, & + ispec,doubling_index, & + xstore(:,:,:,ispec),ystore(:,:,:,ispec),zstore(:,:,:,ispec), & + iboun,iMPIcut_xi,iMPIcut_eta,NPROC_XI,NPROC_ETA, & + UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,NEX_XI,NEX_ETA) + enddo + enddo + enddo + ! end of loop on all the mesh points in current subregion + + else - ! Regular subregion case + ! Irregular subregion case - do ir = ir1,ir2,dir - do iy = iy1,iy2,diy - do ix = ix1,ix2,dix - ! loop over the NGNOD_EIGHT_CORNERS nodes + do ir = ir1,ir2,dir + do iy = iy1,iy2,diy + do ix = ix1,ix2,dix + ! loop on all the elements in the mesh doubling superbrick + do ispec_superbrick = 1,nspec_sb + ! loop on all the corner nodes of this element do ia = 1,NGNOD_EIGHT_CORNERS ! define topological coordinates of this mesh point - ioffset_x = ix+iax*iaddx(ia) - ioffset_y = iy+iay*iaddy(ia) - ioffset_z = ir+iar*iaddz(ia) + ioffset_x = int(ix + iax*x_superbrick(ibool_superbrick(ia,ispec_superbrick))) + ioffset_y = int(iy + iay*y_superbrick(ibool_superbrick(ia,ispec_superbrick))) + ioffset_z = int(ir + iar*z_superbrick(ibool_superbrick(ia,ispec_superbrick))) xelm(ia) = xgrid(ioffset_z,ioffset_x,ioffset_y) yelm(ia) = ygrid(ioffset_z,ioffset_x,ioffset_y) @@ -460,15 +493,12 @@ subroutine cmm_create_mesh_elements() call exit_MPI(myrank,'ispec greater than nspec in mesh creation') endif - ! check if element is on topography - if ((ir == ir2) .and. (isubregion == nmeshregions)) then - ! sets index to topography layer - doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY - else - ! sets index to deeper, basement layers - doubling_index = IFLAG_BASEMENT_TOPO - endif - + ! sets index to deeper, basement layers + ! note: doubling layer can not be at surface, + ! since modulo 2 is zero for doubling layer, it will always have one layer below and one layer above + ! for example: 1 doubling layer -> mesh regions = 3: + ! region 1 is bottom, region 2 is doubling, region 3 top + doubling_index = IFLAG_BASEMENT_TOPO ispec_material_id(ispec) = material_num(ir,ix,iy) ! store coordinates @@ -480,323 +510,303 @@ subroutine cmm_create_mesh_elements() xstore(:,:,:,ispec),ystore(:,:,:,ispec),zstore(:,:,:,ispec), & iboun,iMPIcut_xi,iMPIcut_eta,NPROC_XI,NPROC_ETA, & UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,NEX_XI,NEX_ETA) + enddo enddo enddo - ! end of loop on all the mesh points in current subregion + enddo + ! end of loop on all the mesh points in current subregion + + endif ! regular/doubling region + + ! end of loop on all the subregions of the current region the mesh + enddo + + ! mesh dimensions + x_min = minval(xstore(:,:,:,:)) + x_max = maxval(xstore(:,:,:,:)) + y_min = minval(ystore(:,:,:,:)) + y_max = maxval(ystore(:,:,:,:)) + z_min = minval(zstore(:,:,:,:)) + z_max = maxval(zstore(:,:,:,:)) + call min_all_dp(x_min,x_min_glob) + call max_all_dp(x_max,x_max_glob) + call min_all_dp(y_min,y_min_glob) + call max_all_dp(y_max,y_max_glob) + call min_all_dp(z_min,z_min_glob) + call max_all_dp(z_max,z_max_glob) + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) 'mesh dimensions:' + write(IMAIN,*) ' Xmin and Xmax of the model = ',sngl(x_min_glob),sngl(x_max_glob) + write(IMAIN,*) ' Ymin and Ymax of the model = ',sngl(y_min_glob),sngl(y_max_glob) + write(IMAIN,*) ' Zmin and Zmax of the model = ',sngl(z_min_glob),sngl(z_max_glob) + endif - else + ! compare to exact theoretical value (bottom is always flat) + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) 'exact area = ',sngl((UTM_Y_MAX-UTM_Y_MIN)*(UTM_X_MAX-UTM_X_MIN)),'(m^2)' + write(IMAIN,*) ' = ',sngl((UTM_Y_MAX-UTM_Y_MIN)*(UTM_X_MAX-UTM_X_MIN)/1000./1000.),'(km^2)' + call flush_IMAIN() + endif - ! Irregular subregion case - - do ir = ir1,ir2,dir - do iy = iy1,iy2,diy - do ix = ix1,ix2,dix - ! loop on all the elements in the mesh doubling superbrick - do ispec_superbrick = 1,nspec_sb - ! loop on all the corner nodes of this element - do ia = 1,NGNOD_EIGHT_CORNERS - ! define topological coordinates of this mesh point - ioffset_x = int(ix + iax*x_superbrick(ibool_superbrick(ia,ispec_superbrick))) - ioffset_y = int(iy + iay*y_superbrick(ibool_superbrick(ia,ispec_superbrick))) - ioffset_z = int(ir + iar*z_superbrick(ibool_superbrick(ia,ispec_superbrick))) - - xelm(ia) = xgrid(ioffset_z,ioffset_x,ioffset_y) - yelm(ia) = ygrid(ioffset_z,ioffset_x,ioffset_y) - zelm(ia) = zgrid(ioffset_z,ioffset_x,ioffset_y) - enddo - - ! add one spectral element to the list and store its material number - ispec = ispec + 1 - - if (ispec > nspec) then - call exit_MPI(myrank,'ispec greater than nspec in mesh creation') - endif - - ! sets index to deeper, basement layers - ! note: doubling layer can not be at surface, - ! since modulo 2 is zero for doubling layer, it will always have one layer below and one layer above - ! for example: 1 doubling layer -> mesh regions = 3: - ! region 1 is bottom, region 2 is doubling, region 3 top - doubling_index = IFLAG_BASEMENT_TOPO - ispec_material_id(ispec) = material_num(ir,ix,iy) - - ! store coordinates - call store_coords(xstore,ystore,zstore,xelm,yelm,zelm,ispec,nspec) - - ! detect mesh boundaries - call get_flags_boundaries(nspec,iproc_xi_current,iproc_eta_current, & - ispec,doubling_index, & - xstore(:,:,:,ispec),ystore(:,:,:,ispec),zstore(:,:,:,ispec), & - iboun,iMPIcut_xi,iMPIcut_eta,NPROC_XI,NPROC_ETA, & - UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,NEX_XI,NEX_ETA) + ! check total number of spectral elements created + if (ispec /= nspec) call exit_MPI(myrank,'ispec should equal nspec') + + ! check that we assign a material to each element + if (any(ispec_material_id(:) == -1000)) stop 'Element of undefined material found' + + ! element size + elemsize_min_glob = HUGEVAL + elemsize_max_glob = -HUGEVAL + do ispec = 1,nspec + elemsize_min = HUGEVAL + elemsize_max = -HUGEVAL + ! loops over the four edges that are along X + i1 = 1 + i2 = NGLLX_M + do k = 1, NGLLZ_M, NGLLZ_M-1 + do j = 1, NGLLY_M, NGLLY_M-1 + x1 = xstore(i1,j,k,ispec) + y1 = ystore(i1,j,k,ispec) + z1 = zstore(i1,j,k,ispec) - enddo - enddo - enddo - enddo - ! end of loop on all the mesh points in current subregion + x2 = xstore(i2,j,k,ispec) + y2 = ystore(i2,j,k,ispec) + z2 = zstore(i2,j,k,ispec) - endif ! regular/doubling region + dist = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2) - ! end of loop on all the subregions of the current region the mesh + if (dist < elemsize_min) elemsize_min = dist + if (dist > elemsize_max) elemsize_max = dist + enddo + enddo + ! loops over the four edges that are along Y + j1 = 1 + j2 = NGLLY_M + do k = 1, NGLLZ_M, NGLLZ_M-1 + do i = 1, NGLLX_M, NGLLX_M-1 + x1 = xstore(i,j1,k,ispec) + y1 = ystore(i,j1,k,ispec) + z1 = zstore(i,j1,k,ispec) + + x2 = xstore(i,j2,k,ispec) + y2 = ystore(i,j2,k,ispec) + z2 = zstore(i,j2,k,ispec) + + dist = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2) + + if (dist < elemsize_min) elemsize_min = dist + if (dist > elemsize_max) elemsize_max = dist + enddo + enddo + ! loops over the four edges that are along Z + k1 = 1 + k2 = NGLLZ_M + do j = 1, NGLLY_M, NGLLY_M-1 + do i = 1, NGLLX_M, NGLLX_M-1 + x1 = xstore(i,j,k1,ispec) + y1 = ystore(i,j,k1,ispec) + z1 = zstore(i,j,k1,ispec) + + x2 = xstore(i,j,k2,ispec) + y2 = ystore(i,j,k2,ispec) + z2 = zstore(i,j,k2,ispec) + + dist = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2) + + if (dist < elemsize_min) elemsize_min = dist + if (dist > elemsize_max) elemsize_max = dist + enddo enddo + elemsize_min = sqrt( elemsize_min ) + elemsize_max = sqrt( elemsize_max ) + + elemsize_min_glob = min(elemsize_min_glob, elemsize_min) + elemsize_max_glob = max(elemsize_max_glob, elemsize_max) + enddo + ! element size + elemsize_min = elemsize_min_glob + elemsize_max = elemsize_max_glob + call min_all_cr(elemsize_min,elemsize_min_glob) + call max_all_cr(elemsize_max,elemsize_max_glob) + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) ' Max element size = ',elemsize_max_glob,'(m)' + write(IMAIN,*) ' Min element size = ',elemsize_min_glob,'(m)' + write(IMAIN,*) ' Max/min ratio = ',elemsize_max_glob/elemsize_min_glob + write(IMAIN,*) + call flush_IMAIN() + endif - ! mesh dimensions - x_min = minval(xstore(:,:,:,:)) - x_max = maxval(xstore(:,:,:,:)) - y_min = minval(ystore(:,:,:,:)) - y_max = maxval(ystore(:,:,:,:)) - z_min = minval(zstore(:,:,:,:)) - z_max = maxval(zstore(:,:,:,:)) - call min_all_dp(x_min,x_min_glob) - call max_all_dp(x_max,x_max_glob) - call min_all_dp(y_min,y_min_glob) - call max_all_dp(y_max,y_max_glob) - call min_all_dp(z_min,z_min_glob) - call max_all_dp(z_max,z_max_glob) - if (myrank == 0) then - write(IMAIN,*) - write(IMAIN,*) 'mesh dimensions:' - write(IMAIN,*) ' Xmin and Xmax of the model = ',sngl(x_min_glob),sngl(x_max_glob) - write(IMAIN,*) ' Ymin and Ymax of the model = ',sngl(y_min_glob),sngl(y_max_glob) - write(IMAIN,*) ' Zmin and Zmax of the model = ',sngl(z_min_glob),sngl(z_max_glob) - endif + deallocate(material_num) - ! compare to exact theoretical value (bottom is always flat) - if (myrank == 0) then - write(IMAIN,*) - write(IMAIN,*) 'exact area = ',sngl((UTM_Y_MAX-UTM_Y_MIN)*(UTM_X_MAX-UTM_X_MIN)),'(m^2)' - write(IMAIN,*) ' = ',sngl((UTM_Y_MAX-UTM_Y_MIN)*(UTM_X_MAX-UTM_X_MIN)/1000./1000.),'(km^2)' - call flush_IMAIN() - endif + end subroutine cmm_create_mesh_elements - ! check total number of spectral elements created - if (ispec /= nspec) call exit_MPI(myrank,'ispec should equal nspec') +! +!---------------------------------------------------------------------------------------------------- +! - ! check that we assign a material to each element - if (any(ispec_material_id(:) == -1000)) stop 'Element of undefined material found' + subroutine cmm_create_addressing(nglob) - ! element size - elemsize_min_glob = HUGEVAL - elemsize_max_glob = -HUGEVAL - do ispec = 1,nspec - elemsize_min = HUGEVAL - elemsize_max = -HUGEVAL - ! loops over the four edges that are along X - i1 = 1 - i2 = NGLLX_M - do k = 1, NGLLZ_M, NGLLZ_M-1 - do j = 1, NGLLY_M, NGLLY_M-1 - x1 = xstore(i1,j,k,ispec) - y1 = ystore(i1,j,k,ispec) - z1 = zstore(i1,j,k,ispec) - - x2 = xstore(i2,j,k,ispec) - y2 = ystore(i2,j,k,ispec) - z2 = zstore(i2,j,k,ispec) - - dist = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2) - - if (dist < elemsize_min) elemsize_min = dist - if (dist > elemsize_max) elemsize_max = dist - enddo - enddo - ! loops over the four edges that are along Y - j1 = 1 - j2 = NGLLY_M - do k = 1, NGLLZ_M, NGLLZ_M-1 - do i = 1, NGLLX_M, NGLLX_M-1 - x1 = xstore(i,j1,k,ispec) - y1 = ystore(i,j1,k,ispec) - z1 = zstore(i,j1,k,ispec) - - x2 = xstore(i,j2,k,ispec) - y2 = ystore(i,j2,k,ispec) - z2 = zstore(i,j2,k,ispec) - - dist = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2) - - if (dist < elemsize_min) elemsize_min = dist - if (dist > elemsize_max) elemsize_max = dist - enddo - enddo - ! loops over the four edges that are along Z - k1 = 1 - k2 = NGLLZ_M - do j = 1, NGLLY_M, NGLLY_M-1 - do i = 1, NGLLX_M, NGLLX_M-1 - x1 = xstore(i,j,k1,ispec) - y1 = ystore(i,j,k1,ispec) - z1 = zstore(i,j,k1,ispec) + use constants, only: NDIM,IOVTK,IMAIN + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NGLLCUBE_M - x2 = xstore(i,j,k2,ispec) - y2 = ystore(i,j,k2,ispec) - z2 = zstore(i,j,k2,ispec) + use meshfem3D_par, only: myrank,prname,ibool,xstore,ystore,zstore, & + UTM_X_MIN,UTM_X_MAX,NGLOB_AB,nspec,npointot - dist = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2) + use create_meshfem_par - if (dist < elemsize_min) elemsize_min = dist - if (dist > elemsize_max) elemsize_max = dist - enddo - enddo - elemsize_min = sqrt( elemsize_min ) - elemsize_max = sqrt( elemsize_max ) + implicit none - elemsize_min_glob = min(elemsize_min_glob, elemsize_min) - elemsize_max_glob = max(elemsize_max_glob, elemsize_max) - enddo - ! element size - elemsize_min = elemsize_min_glob - elemsize_max = elemsize_max_glob - call min_all_cr(elemsize_min,elemsize_min_glob) - call max_all_cr(elemsize_max,elemsize_max_glob) - ! user output - if (myrank == 0) then - write(IMAIN,*) - write(IMAIN,*) ' Max element size = ',elemsize_max_glob,'(m)' - write(IMAIN,*) ' Min element size = ',elemsize_min_glob,'(m)' - write(IMAIN,*) ' Max/min ratio = ',elemsize_max_glob/elemsize_min_glob - write(IMAIN,*) - call flush_IMAIN() - endif + integer, intent(out) :: nglob - end subroutine cmm_create_mesh_elements + ! local parameters + integer :: i,j,k,ispec,ier + integer :: ieoff,ilocnum -! -!---------------------------------------------------------------------------------------------------- -! + ! variables for creating array ibool (some arrays also used for AVS or DX files) + integer, dimension(:), allocatable :: iglob,locval + logical, dimension(:), allocatable :: ifseg + double precision, dimension(:), allocatable :: xp,yp,zp - subroutine cmm_create_addressing() + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) 'creating indirect addressing for unstructured mesh' + write(IMAIN,*) + call flush_IMAIN() + endif - use meshfem3D_par, only: prname + ! allocate memory for arrays + allocate(iglob(npointot),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1272') + allocate(locval(npointot),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1273') + allocate(ifseg(npointot),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1274') + allocate(xp(npointot),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1275') + allocate(yp(npointot),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1276') + allocate(zp(npointot),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'error allocating array 1277') + + ! puts x,y,z locations into 1D arrays + do ispec = 1,nspec + ieoff = NGLLCUBE_M*(ispec-1) + ilocnum = 0 + do k = 1,NGLLZ_M + do j = 1,NGLLY_M + do i = 1,NGLLX_M + ilocnum = ilocnum + 1 + xp(ilocnum+ieoff) = xstore(i,j,k,ispec) + yp(ilocnum+ieoff) = ystore(i,j,k,ispec) + zp(ilocnum+ieoff) = zstore(i,j,k,ispec) + enddo + enddo + enddo + enddo - implicit none - ! local parameters - integer :: i,j,k,ispec,ier - integer :: ieoff,ilocnum + ! sorts xp,yp,zp in lexicographical order (increasing values) + call get_global(npointot,xp,yp,zp,iglob,locval,ifseg,nglob,UTM_X_MIN,UTM_X_MAX) + ! checks nglob range with pre-computed values + ! note: if mesh squeezes elements such that we can't distinguish two close-by mesh points anymore + ! the total number of mesh points might have changed + if (nglob /= NGLOB_AB) then ! user output + print *,'Error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB if (myrank == 0) then write(IMAIN,*) - write(IMAIN,*) 'creating indirect addressing for unstructured mesh' + write(IMAIN,*) 'Error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB write(IMAIN,*) + write(IMAIN,*) 'writing out problematic mesh: mesh_debug.vtk' + write(IMAIN,*) 'please check your mesh setup...' call flush_IMAIN() endif - ! puts x,y,z locations into 1D arrays + ! debug file output + ! vtk file format output + open(IOVTK,file=prname(1:len_trim(prname))//'mesh_debug.vtk',status='unknown') + write(IOVTK,'(a)') '# vtk DataFile Version 3.1' + write(IOVTK,'(a)') 'material model VTK file' + write(IOVTK,'(a)') 'ASCII' + write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID' + write(IOVTK, '(a,i12,a)') 'POINTS ', nspec*NGLLX_M*NGLLY_M*NGLLZ_M, ' float' + ilocnum = 0 + ibool(:,:,:,:) = 0 do ispec = 1,nspec - ieoff = NGLLCUBE_M*(ispec-1) - ilocnum = 0 do k = 1,NGLLZ_M do j = 1,NGLLY_M do i = 1,NGLLX_M ilocnum = ilocnum + 1 - xp(ilocnum+ieoff) = xstore(i,j,k,ispec) - yp(ilocnum+ieoff) = ystore(i,j,k,ispec) - zp(ilocnum+ieoff) = zstore(i,j,k,ispec) + ibool(i,j,k,ispec) = ilocnum + write(IOVTK,*) xstore(i,j,k,ispec),ystore(i,j,k,ispec),zstore(i,j,k,ispec) enddo enddo enddo enddo - - ! sorts xp,yp,zp in lexicographical order (increasing values) - call get_global(npointot,xp,yp,zp,iglob,locval,ifseg,nglob,UTM_X_MIN,UTM_X_MAX) - - ! checks nglob range with pre-computed values - ! note: if mesh squeezes elements such that we can't distinguish two close-by mesh points anymore - ! the total number of mesh points might have changed - if (nglob /= NGLOB_AB) then - ! user output - print *,'Error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB - if (myrank == 0) then - write(IMAIN,*) - write(IMAIN,*) 'Error nglob: sorted value ',nglob,'differs from pre-computed ',NGLOB_AB - write(IMAIN,*) - write(IMAIN,*) 'writing out problematic mesh: mesh_debug.vtk' - write(IMAIN,*) 'please check your mesh setup...' - call flush_IMAIN() - endif - - ! debug file output - ! vtk file format output - open(IOVTK,file=prname(1:len_trim(prname))//'mesh_debug.vtk',status='unknown') - write(IOVTK,'(a)') '# vtk DataFile Version 3.1' - write(IOVTK,'(a)') 'material model VTK file' - write(IOVTK,'(a)') 'ASCII' - write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID' - write(IOVTK, '(a,i12,a)') 'POINTS ', nspec*NGLLX_M*NGLLY_M*NGLLZ_M, ' float' - ilocnum = 0 - ibool(:,:,:,:) = 0 - do ispec = 1,nspec - do k = 1,NGLLZ_M - do j = 1,NGLLY_M - do i = 1,NGLLX_M - ilocnum = ilocnum + 1 - ibool(i,j,k,ispec) = ilocnum - write(IOVTK,*) xstore(i,j,k,ispec),ystore(i,j,k,ispec),zstore(i,j,k,ispec) - enddo - enddo - enddo - enddo - write(IOVTK,*) - ! note: indices for vtk start at 0 - write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9 - do ispec = 1,nspec - write(IOVTK,'(9i12)') 8, & - ibool(1,1,1,ispec)-1,ibool(2,1,1,ispec)-1,ibool(2,2,1,ispec)-1,ibool(1,2,1,ispec)-1, & - ibool(1,1,2,ispec)-1,ibool(2,1,2,ispec)-1,ibool(2,2,2,ispec)-1,ibool(1,2,2,ispec)-1 - enddo - ibool(:,:,:,:) = 0 - write(IOVTK,*) - ! type: hexahedrons - write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec - write(IOVTK,'(6i12)') (12,ispec=1,nspec) - write(IOVTK,*) - write(IOVTK,'(a,i12)') "CELL_DATA ",nspec - write(IOVTK,'(a)') "SCALARS elem_val float" - write(IOVTK,'(a)') "LOOKUP_TABLE default" - do ispec = 1,nspec - write(IOVTK,*) ispec_material_id(ispec) - enddo - write(IOVTK,*) '' - close(IOVTK) - - ! stop mesher - call exit_MPI(myrank,'incorrect global number, please check mesh input parameters') - endif - - ! put in classical format - allocate(nodes_coords(nglob,NDIM),stat=ier) - if (ier /= 0) call exit_MPI_without_rank('error allocating array 1280') - if (ier /= 0) stop 'Error allocating array nodes_coords' - - nodes_coords(:,:) = 0.0d0 + write(IOVTK,*) + ! note: indices for vtk start at 0 + write(IOVTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9 + do ispec = 1,nspec + write(IOVTK,'(9i12)') 8, & + ibool(1,1,1,ispec)-1,ibool(2,1,1,ispec)-1,ibool(2,2,1,ispec)-1,ibool(1,2,1,ispec)-1, & + ibool(1,1,2,ispec)-1,ibool(2,1,2,ispec)-1,ibool(2,2,2,ispec)-1,ibool(1,2,2,ispec)-1 + enddo ibool(:,:,:,:) = 0 - + write(IOVTK,*) + ! type: hexahedrons + write(IOVTK,'(a,i12)') "CELL_TYPES ",nspec + write(IOVTK,'(6i12)') (12,ispec=1,nspec) + write(IOVTK,*) + write(IOVTK,'(a,i12)') "CELL_DATA ",nspec + write(IOVTK,'(a)') "SCALARS elem_val float" + write(IOVTK,'(a)') "LOOKUP_TABLE default" do ispec = 1,nspec - ieoff = NGLLCUBE_M*(ispec-1) - ilocnum = 0 - do k = 1,NGLLZ_M - do j = 1,NGLLY_M - do i = 1,NGLLX_M - ilocnum = ilocnum + 1 - ibool(i,j,k,ispec) = iglob(ilocnum+ieoff) - nodes_coords(iglob(ilocnum+ieoff),1) = xstore(i,j,k,ispec) - nodes_coords(iglob(ilocnum+ieoff),2) = ystore(i,j,k,ispec) - nodes_coords(iglob(ilocnum+ieoff),3) = zstore(i,j,k,ispec) - enddo + write(IOVTK,*) ispec_material_id(ispec) + enddo + write(IOVTK,*) '' + close(IOVTK) + + ! stop mesher + call exit_MPI(myrank,'incorrect global number, please check mesh input parameters') + endif + + ! put in classical format + allocate(nodes_coords(nglob,NDIM),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 1280') + if (ier /= 0) stop 'Error allocating array nodes_coords' + + nodes_coords(:,:) = 0.0d0 + ibool(:,:,:,:) = 0 + + do ispec = 1,nspec + ieoff = NGLLCUBE_M*(ispec-1) + ilocnum = 0 + do k = 1,NGLLZ_M + do j = 1,NGLLY_M + do i = 1,NGLLX_M + ilocnum = ilocnum + 1 + ibool(i,j,k,ispec) = iglob(ilocnum+ieoff) + nodes_coords(iglob(ilocnum+ieoff),1) = xstore(i,j,k,ispec) + nodes_coords(iglob(ilocnum+ieoff),2) = ystore(i,j,k,ispec) + nodes_coords(iglob(ilocnum+ieoff),3) = zstore(i,j,k,ispec) enddo enddo enddo + enddo - ! checks ibool range - if (minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= nglob) then - print *,'Error ibool: maximum value ',maxval(ibool(:,:,:,:)) ,'should be ',nglob - call exit_MPI(myrank,'incorrect global ibool numbering') - endif + ! checks ibool range + if (minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= nglob) then + print *,'Error ibool: maximum value ',maxval(ibool(:,:,:,:)) ,'should be ',nglob + call exit_MPI(myrank,'incorrect global ibool numbering') + endif - end subroutine cmm_create_addressing + deallocate(iglob,locval,ifseg,xp,yp,zp) + end subroutine cmm_create_addressing - end subroutine create_meshfem_mesh diff --git a/src/meshfem3D/create_visual_files.f90 b/src/meshfem3D/create_visual_files.f90 index d437fa0ce..67a152d3b 100644 --- a/src/meshfem3D/create_visual_files.f90 +++ b/src/meshfem3D/create_visual_files.f90 @@ -30,11 +30,10 @@ subroutine create_visual_files(CREATE_ABAQUS_FILES,CREATE_DX_FILES,CREATE_VTK_FI nspec,nglob,prname,nodes_coords,ibool,ispec_material_id) use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M implicit none - include "constants_meshfem3D.h" - ! Mesh files for visualization logical :: CREATE_ABAQUS_FILES,CREATE_DX_FILES,CREATE_VTK_FILES diff --git a/src/meshfem3D/define_subregions_heuristic.f90 b/src/meshfem3D/define_subregions_heuristic.f90 index ce7976491..5452366b9 100644 --- a/src/meshfem3D/define_subregions_heuristic.f90 +++ b/src/meshfem3D/define_subregions_heuristic.f90 @@ -37,11 +37,10 @@ subroutine define_subregions_heuristic(myrank,isubregion,iaddx,iaddy,iaddz, & ! to 120 degrees in doubling regions use constants + use constants_meshfem3D, only: ITYPE_UNUSUAL_1,ITYPE_UNUSUAL_1p,ITYPE_UNUSUAL_4,ITYPE_UNUSUAL_4p implicit none - include "constants_meshfem3D.h" - integer myrank integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir integer iax,iay,iar diff --git a/src/meshfem3D/define_superbrick.f90 b/src/meshfem3D/define_superbrick.f90 index ab2452083..48636243f 100644 --- a/src/meshfem3D/define_superbrick.f90 +++ b/src/meshfem3D/define_superbrick.f90 @@ -31,11 +31,10 @@ subroutine define_superbrick(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb) use constants + use constants_meshfem3D, only: NGLOB_DOUBLING_SUPERBRICK,NSPEC_DOUBLING_SUPERBRICK implicit none - include "constants_meshfem3D.h" - integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick logical, dimension(NSPEC_DOUBLING_SUPERBRICK,6) :: iboun_sb @@ -661,11 +660,10 @@ end subroutine define_superbrick subroutine define_superbrick_one_layer(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb) use constants + use constants_meshfem3D, only: NGLOB_DOUBLING_SUPERBRICK,NSPEC_DOUBLING_SUPERBRICK implicit none - include "constants_meshfem3D.h" - integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick logical, dimension(NSPEC_DOUBLING_SUPERBRICK,6) :: iboun_sb @@ -1209,11 +1207,10 @@ end subroutine define_superbrick_one_layer subroutine define_basic_doubling_brick(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb,case_num) use constants + use constants_meshfem3D, only: NGLOB_DOUBLING_SUPERBRICK,NSPEC_DOUBLING_SUPERBRICK implicit none - include "constants_meshfem3D.h" - integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick logical, dimension(NSPEC_DOUBLING_SUPERBRICK,6) :: iboun_sb diff --git a/src/meshfem3D/determine_cavity.f90 b/src/meshfem3D/determine_cavity.f90 index 5696aa7e5..04e17060d 100644 --- a/src/meshfem3D/determine_cavity.f90 +++ b/src/meshfem3D/determine_cavity.f90 @@ -28,15 +28,13 @@ subroutine cmm_determine_cavity(nglob) + use constants, only: MF_IN_DATA_FILES,MAX_STRING_LEN,IMAIN,HUGEVAL,TINYVAL,NDIM + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M + use create_meshfem_par, only: nodes_coords,ispec_material_id,iboun,iMPIcut_xi,iMPIcut_eta use meshfem3D_par, only: ibool,xstore,ystore,zstore,nspec,NPROC_XI,NPROC_ETA,myrank,CAVITY_FILE - ! create the different regions of the mesh - use constants, only: MF_IN_DATA_FILES,MAX_STRING_LEN,IMAIN,CUSTOM_REAL,HUGEVAL,TINYVAL,NDIM - - use constants_meshfem3D, only: NGLLCUBE_M,NGLLX_M,NGLLY_M,NGLLZ_M - implicit none integer,intent(inout) :: nglob @@ -70,10 +68,8 @@ subroutine cmm_determine_cavity(nglob) character(len=MAX_STRING_LEN) :: filename - !-------------------------------cavity---------------------------------------- - ! begin cavity ! default ncavity = 0 diff --git a/src/meshfem3D/earth_chunk.f90 b/src/meshfem3D/earth_chunk.f90 index 2dfdc5ef4..58f2feae9 100644 --- a/src/meshfem3D/earth_chunk.f90 +++ b/src/meshfem3D/earth_chunk.f90 @@ -31,7 +31,7 @@ subroutine earth_chunk_HEX8_Mesher(NGNOD) - use constants, only: NGLLX, NGLLY, NGLLZ, NDIM, R_EARTH, PI, ZERO, TINYVAL, & + use constants, only: NGLLX, NGLLY, NGLLZ, NDIM, R_EARTH, ZERO, & old_DSM_coupling_from_Vadim, INJECTION_TECHNIQUE_IS_AXISEM, INJECTION_TECHNIQUE_IS_DSM use shared_parameters, only: INJECTION_TECHNIQUE_TYPE @@ -757,7 +757,7 @@ end subroutine earth_chunk_HEX8_Mesher subroutine earth_chunk_HEX27_Mesher(NGNOD) - use constants, only: NGLLX, NGLLY, NGLLZ, NDIM, R_EARTH, PI, ZERO, TINYVAL, & + use constants, only: NGLLX, NGLLY, NGLLZ, NDIM, R_EARTH, ZERO, & old_DSM_coupling_from_Vadim, INJECTION_TECHNIQUE_IS_AXISEM, INJECTION_TECHNIQUE_IS_DSM use shared_parameters, only: INJECTION_TECHNIQUE_TYPE diff --git a/src/meshfem3D/get_MPI_cutplanes_eta.f90 b/src/meshfem3D/get_MPI_cutplanes_eta.f90 index a2cd662c4..ac8105d16 100644 --- a/src/meshfem3D/get_MPI_cutplanes_eta.f90 +++ b/src/meshfem3D/get_MPI_cutplanes_eta.f90 @@ -38,11 +38,10 @@ subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, & ! in the solver except if we want to have periodic conditions use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M implicit none - include "constants_meshfem3D.h" - integer nspec,myrank integer NSPEC2D_A_XI,NSPEC2D_B_XI diff --git a/src/meshfem3D/get_MPI_cutplanes_xi.f90 b/src/meshfem3D/get_MPI_cutplanes_xi.f90 index 713d75467..8593f67df 100644 --- a/src/meshfem3D/get_MPI_cutplanes_xi.f90 +++ b/src/meshfem3D/get_MPI_cutplanes_xi.f90 @@ -38,11 +38,10 @@ subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, & ! in the solver except if we want to have periodic conditions use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M implicit none - include "constants_meshfem3D.h" - integer nspec,myrank integer NSPEC2D_A_ETA,NSPEC2D_B_ETA diff --git a/src/meshfem3D/get_flags_boundaries.f90 b/src/meshfem3D/get_flags_boundaries.f90 index 036305b2d..a5a357cdb 100644 --- a/src/meshfem3D/get_flags_boundaries.f90 +++ b/src/meshfem3D/get_flags_boundaries.f90 @@ -31,11 +31,10 @@ subroutine get_flags_boundaries(nspec,iproc_xi,iproc_eta,ispec,idoubling, & UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,NEX_XI,NEX_ETA) use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,IFLAG_ONE_LAYER_TOPOGRAPHY implicit none - include "constants_meshfem3D.h" - integer,intent(in) :: nspec integer,intent(in) :: ispec,idoubling integer,intent(in) :: NPROC_XI,NPROC_ETA diff --git a/src/meshfem3D/meshfem3D_adios_stubs.f90 b/src/meshfem3D/meshfem3D_adios_stubs.f90 index f4162ad20..a6c38cba8 100644 --- a/src/meshfem3D/meshfem3D_adios_stubs.f90 +++ b/src/meshfem3D/meshfem3D_adios_stubs.f90 @@ -45,12 +45,12 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & nspec_CPML,CPML_to_spec,CPML_regions,is_CPML) use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES + use adios_manager_mod implicit none - include "constants_meshfem3D.h" - character(len=MAX_STRING_LEN) :: LOCAL_PATH integer :: myrank, sizeprocs integer :: nspec diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk index e894d9f0b..398b4cc82 100644 --- a/src/meshfem3D/rules.mk +++ b/src/meshfem3D/rules.mk @@ -58,7 +58,7 @@ meshfem3D_OBJECTS = \ $O/get_MPI_cutplanes_eta.mesh.o \ $O/get_MPI_cutplanes_xi.mesh.o \ $O/meshfem3D.mesh.o \ - $O/meshfem3D_par.mesh.o \ + $O/meshfem3D_par.mesh_module.o \ $O/read_mesh_parameter_file.mesh.o \ $O/read_value_mesh_parameters.mesh.o \ $O/save_databases.mesh.o \ @@ -150,12 +150,7 @@ $E/xmeshfem3D: $(XMESHFEM_OBJECTS) ### Module dependencies ### - -$O/meshfem3D.mesh.o: $O/meshfem3D_par.mesh.o $O/chunk_earth_mesh_mod.mesh.o -$O/create_meshfem_mesh.mesh.o: $O/meshfem3D_par.mesh.o -$O/create_CPML_regions.mesh.o: $O/meshfem3D_par.mesh.o -$O/create_interfaces_mesh.mesh.o: $O/meshfem3D_par.mesh.o -$O/read_mesh_parameter_file.mesh.o: $O/meshfem3D_par.mesh.o +$O/meshfem3D.mesh.o: $O/chunk_earth_mesh_mod.mesh.o $O/determine_cavity.mesh.o: $O/create_meshfem_mesh.mesh.o ## adios @@ -173,10 +168,13 @@ $O/adios_helpers.shared_adios.o: \ #### rule to build each .o file below #### -$O/%.mesh.o: $S/%.f90 $O/shared_par.shared_module.o +$O/%.mesh_module.o: $S/%.f90 $O/shared_par.shared_module.o + ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< + +$O/%.mesh.o: $S/%.f90 $O/shared_par.shared_module.o $O/meshfem3D_par.mesh_module.o ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< -$O/%.mesh.o: $S/%.F90 $O/shared_par.shared_module.o +$O/%.mesh.o: $S/%.F90 $O/shared_par.shared_module.o $O/meshfem3D_par.mesh_module.o ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< @@ -184,10 +182,10 @@ $O/%.mesh.o: $S/%.F90 $O/shared_par.shared_module.o ### ADIOS compilation ### -$O/%.mesh_adios.o: $S/%.F90 $O/shared_par.shared_module.o +$O/%.mesh_adios.o: $S/%.F90 $O/shared_par.shared_module.o $O/meshfem3D_par.mesh_module.o ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< -$O/%.mesh_adios.o: $S/%.f90 $O/shared_par.shared_module.o +$O/%.mesh_adios.o: $S/%.f90 $O/shared_par.shared_module.o $O/meshfem3D_par.mesh_module.o ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< $O/%.mesh_noadios.o: $S/%.F90 diff --git a/src/meshfem3D/save_databases.F90 b/src/meshfem3D/save_databases.F90 index f79ea768e..cceb4f38c 100644 --- a/src/meshfem3D/save_databases.F90 +++ b/src/meshfem3D/save_databases.F90 @@ -36,12 +36,12 @@ subroutine save_databases(prname,nspec,nglob,iproc_xi,iproc_eta, & xstore, ystore, zstore) use constants, only: MAX_STRING_LEN,IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,SAVE_MESH_AS_CUBIT,NDIM + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES + use shared_parameters, only: COUPLE_WITH_INJECTION_TECHNIQUE implicit none - include "constants_meshfem3D.h" - ! number of spectral elements in each block integer nspec @@ -426,12 +426,11 @@ subroutine save_output_mesh_files_as_cubit(nspec,nglob, ibool,nodes_coords, ispe NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NMATERIALS,material_properties, & ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top) - use constants, only: MAX_STRING_LEN,IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC, NGLLX, NGLLY, NGLLZ, NDIM, ZERO + use constants, only: NDIM + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES implicit none - include "constants_meshfem3D.h" - integer, parameter :: IIN_database = 15 ! number of spectral elements in each block integer :: nspec @@ -441,7 +440,7 @@ subroutine save_output_mesh_files_as_cubit(nspec,nglob, ibool,nodes_coords, ispe ! MPI Cartesian topology ! E for East (= XI_MIN), W for West (= XI_MAX), S for South (= ETA_MIN), N for North (= ETA_MAX) - integer, parameter :: W=1,E=2,S=3,N=4,NW=5,NE=6,SE=7,SW=8 + !integer, parameter :: W=1,E=2,S=3,N=4,NW=5,NE=6,SE=7,SW=8 ! arrays with the mesh integer :: ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) @@ -579,21 +578,20 @@ subroutine save_output_mesh_files_for_coupled_model(nspec, & ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & xgrid,ygrid,zgrid) - use constants, only: MAX_STRING_LEN,IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC, NGLLX, NGLLY, NGLLZ, NDIM, ZERO, & + use constants, only: NGLLX, NGLLY, NGLLZ, NDIM, ZERO, & INJECTION_TECHNIQUE_IS_AXISEM + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M use shared_parameters, only: NGNOD,COUPLE_WITH_INJECTION_TECHNIQUE,INJECTION_TECHNIQUE_TYPE implicit none - include "constants_meshfem3D.h" - ! number of spectral elements in each block integer :: nspec ! MPI Cartesian topology ! E for East (= XI_MIN), W for West (= XI_MAX), S for South (= ETA_MIN), N for North (= ETA_MAX) - integer, parameter :: W=1,E=2,S=3,N=4,NW=5,NE=6,SE=7,SW=8 + !integer, parameter :: W=1,E=2,S=3,N=4,NW=5,NE=6,SE=7,SW=8 !! VM VM add all GLL points for Axisem coupling double precision, dimension(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) :: xgrid, ygrid, zgrid diff --git a/src/meshfem3D/save_databases_adios.F90 b/src/meshfem3D/save_databases_adios.F90 index 377c5aa29..a663eb38b 100644 --- a/src/meshfem3D/save_databases_adios.F90 +++ b/src/meshfem3D/save_databases_adios.F90 @@ -46,6 +46,7 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & use constants, only: MAX_STRING_LEN,IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,ADIOS_TRANSPORT_METHOD, & NGLLX,NGLLY,NGLLZ,NDIM + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES use adios_helpers_mod, only: define_adios_global_array1d,define_adios_scalar, & write_adios_global_1d_array,write_adios_global_string_1d_array, & @@ -55,8 +56,6 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & implicit none - include "constants_meshfem3D.h" - ! MPI variables integer :: myrank, sizeprocs diff --git a/src/meshfem3D/store_boundaries.f90 b/src/meshfem3D/store_boundaries.f90 index 0dc0711cd..27501bd96 100644 --- a/src/meshfem3D/store_boundaries.f90 +++ b/src/meshfem3D/store_boundaries.f90 @@ -35,8 +35,6 @@ subroutine store_boundaries(myrank,iboun,nspec, & implicit none - include "constants_meshfem3D.h" - integer nspec,myrank integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX @@ -139,11 +137,10 @@ subroutine get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, & NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX) use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M implicit none - include "constants_meshfem3D.h" - integer nspec,myrank integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX @@ -373,8 +370,6 @@ subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D,jacobian2 implicit none - include "constants_meshfem3D.h" - ! generic routine that accepts any polynomial degree in each direction integer ispecb,NGLLA,NGLLB,NSPEC2DMAX_AB,myrank @@ -411,7 +406,7 @@ subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D,jacobian2 uny = zxi*xeta-zeta*xxi unz = xxi*yeta-xeta*yxi jacobian = dsqrt(unx**2+uny**2+unz**2) - if (jacobian == ZERO) call exit_MPI(myrank,'2D Jacobian undefined') + if (jacobian <= ZERO) call exit_MPI(myrank,'2D Jacobian undefined') ! normalize normal vector and store surface jacobian ! distinguish if single or double precision for reals diff --git a/src/meshfem3D/store_coords.f90 b/src/meshfem3D/store_coords.f90 index 62ec95c0c..b03c97616 100644 --- a/src/meshfem3D/store_coords.f90 +++ b/src/meshfem3D/store_coords.f90 @@ -28,11 +28,10 @@ subroutine store_coords(xstore,ystore,zstore,xelm,yelm,zelm,ispec,nspec) use constants + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M implicit none - include "constants_meshfem3D.h" - integer,intent(in) :: ispec,nspec double precision, dimension(NGNOD_EIGHT_CORNERS),intent(in) :: xelm,yelm,zelm diff --git a/src/shared/get_attenuation_model.f90 b/src/shared/get_attenuation_model.f90 index 53602df96..a18feed87 100644 --- a/src/shared/get_attenuation_model.f90 +++ b/src/shared/get_attenuation_model.f90 @@ -948,7 +948,7 @@ subroutine attenuation_invert_by_simplex(t2, t1, n, Q_real, tau_s, tau_eps) integer, parameter :: nf = 100 double precision, dimension(nf) :: f - double precision, parameter :: PI = 3.14159265358979d0 + !double precision, parameter :: PI = 3.14159265358979d0 double precision, external :: attenuation_eval ! Values to be passed into the simplex minimization routine diff --git a/src/shared/gll_library.f90 b/src/shared/gll_library.f90 index 73f0b6269..84e9d225e 100644 --- a/src/shared/gll_library.f90 +++ b/src/shared/gll_library.f90 @@ -606,11 +606,13 @@ double precision function pnglj(z,n) ! !------------------------------------------------------------------------ + use constants, only: TINYVAL,ONE + implicit none - include "constants.h" - double precision z - integer n + double precision,intent(in) :: z + integer,intent(in) :: n + double precision, external :: pnleg if (abs(z+1.d0) > TINYVAL) then ! if (z /= -1.d0) @@ -634,11 +636,13 @@ double precision function pndglj(z,n) ! !------------------------------------------------------------------------ + use constants, only: TINYVAL,ONE + implicit none - include "constants.h" - double precision z - integer n + double precision,intent(in) :: z + integer,intent(in) :: n + double precision, external :: pnleg, pndleg if (abs(z+1.d0) > TINYVAL) then ! if (z /= -1.d0) diff --git a/src/shared/parallel.f90 b/src/shared/parallel.f90 index 9b062c2f6..e3b1f1531 100644 --- a/src/shared/parallel.f90 +++ b/src/shared/parallel.f90 @@ -1432,8 +1432,6 @@ subroutine gatherv_all_i(sendbuf, sendcnt, recvbuf, recvcount, recvoffset,recvco implicit none - include "precision.h" - integer :: sendcnt,recvcounttot,NPROC integer, dimension(NPROC) :: recvcount,recvoffset integer, dimension(sendcnt) :: sendbuf diff --git a/src/shared/read_parameter_file.F90 b/src/shared/read_parameter_file.F90 index fcbdb1694..a18a7b9da 100644 --- a/src/shared/read_parameter_file.F90 +++ b/src/shared/read_parameter_file.F90 @@ -1062,7 +1062,7 @@ subroutine get_number_of_sources(sources_filename) ! determines number of sources depending on number of lines in source file ! (only executed by master process) - use constants, only: IIN,IMAIN,IN_DATA_FILES,HUGEVAL,TINYVAL, & + use constants, only: IIN,IN_DATA_FILES,HUGEVAL,TINYVAL, & NLINES_PER_CMTSOLUTION_SOURCE,NLINES_PER_FORCESOLUTION_SOURCE use shared_parameters diff --git a/src/specfem3D/compute_stacey_viscoelastic.F90 b/src/specfem3D/compute_stacey_viscoelastic.F90 index 52a100939..54dd07e6d 100644 --- a/src/specfem3D/compute_stacey_viscoelastic.F90 +++ b/src/specfem3D/compute_stacey_viscoelastic.F90 @@ -46,7 +46,7 @@ subroutine compute_stacey_viscoelastic(NSPEC_AB,NGLOB_AB,accel, & use shared_parameters, only: COUPLE_WITH_INJECTION_TECHNIQUE,INJECTION_TECHNIQUE_TYPE,RECIPROCITY_AND_KH_INTEGRAL ! added by Ping Tong (TP / Tong Ping) for the FK3D calculation - use specfem_par_coupling, only: npt,nbdglb, NTIME_BETWEEN_FFT, & + use specfem_par_coupling, only: npt,nbdglb, & VX_t, VY_t, VZ_t, TX_t, TY_t, TZ_t, NP_RESAMP, & vx_FK,vy_FK,vz_FK,tx_FK,ty_FK,tz_FK diff --git a/src/specfem3D/get_cmt.f90 b/src/specfem3D/get_cmt.f90 index cc6026f0e..a0c854a4e 100644 --- a/src/specfem3D/get_cmt.f90 +++ b/src/specfem3D/get_cmt.f90 @@ -28,7 +28,7 @@ subroutine get_cmt(CMTSOLUTION,tshift_cmt,hdur,lat,long,depth,moment_tensor, & DT,NSOURCES,min_tshift_cmt_original,user_source_time_function) - use constants, only: IIN,IN_DATA_FILES,MAX_STRING_LEN,CUSTOM_REAL + use constants, only: IIN,MAX_STRING_LEN,CUSTOM_REAL use shared_parameters, only: USE_EXTERNAL_SOURCE_FILE,NSTEP_STF,NSOURCES_STF implicit none diff --git a/src/specfem3D/get_force.f90 b/src/specfem3D/get_force.f90 index 3dcfe78b7..e893b63a2 100644 --- a/src/specfem3D/get_force.f90 +++ b/src/specfem3D/get_force.f90 @@ -31,7 +31,7 @@ subroutine get_force(FORCESOLUTION,tshift_force,hdur,lat,long,depth,NSOURCES, & comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP, & user_source_time_function) - use constants, only: IIN,IN_DATA_FILES,MAX_STRING_LEN,TINYVAL,CUSTOM_REAL + use constants, only: IIN,MAX_STRING_LEN,TINYVAL,CUSTOM_REAL use shared_parameters, only: USE_EXTERNAL_SOURCE_FILE,NSTEP_STF,NSOURCES_STF implicit none diff --git a/src/specfem3D/make_gravity.f90 b/src/specfem3D/make_gravity.f90 index 6beb00bf4..a31df2a91 100644 --- a/src/specfem3D/make_gravity.f90 +++ b/src/specfem3D/make_gravity.f90 @@ -458,112 +458,118 @@ subroutine deriv(y,yprime,n,r,ndis,kdis,s1,s2,s3) implicit none ! Argument variables - integer kdis(28),n,ndis - double precision r(n),s1(n),s2(n),s3(n) - double precision y(n),yprime(n) + integer :: kdis(28),n,ndis + double precision :: r(n),s1(n),s2(n),s3(n) + double precision :: y(n),yprime(n) ! Local variables - integer i,j,j1,j2 - integer k,nd,ndp - double precision a0,b0,b1 - double precision f(3,1000),h,h2,h2a - double precision h2b,h3a,ha,s13 - double precision s21,s32,yy(3) + integer :: i,j,j1,j2 + integer :: k,nd,ndp + double precision :: a0,b0,b1 + double precision :: f(3,1000),h,h2,h2a + double precision :: h2b,h3a,ha,s13 + double precision :: s21,s32,yy(3) yy(1) = 0.d0 yy(2) = 0.d0 yy(3) = 0.d0 - ndp=ndis+1 - do 3 nd=1,ndp - if (nd == 1) goto 4 - if (nd == ndp) goto 5 - j1=kdis(nd-1)+1 - j2=kdis(nd)-2 - goto 6 - 4 j1=1 - j2=kdis(1)-2 - goto 6 - 5 j1=kdis(ndis)+1 - j2=n-2 - 6 if ((j2+1-j1) > 0) goto 11 - j2=j2+2 - yy(1)=(y(j2)-y(j1))/(r(j2)-r(j1)) - s1(j1)=yy(1) - s1(j2)=yy(1) - s2(j1)=yy(2) - s2(j2)=yy(2) - s3(j1)=yy(3) - s3(j2)=yy(3) - goto 3 - 11 a0=0.0d0 - if (j1 == 1) goto 7 - h=r(j1+1)-r(j1) - h2=r(j1+2)-r(j1) - yy(1)=h*h2*(h2-h) - h=h*h - h2=h2*h2 - b0=(y(j1)*(h-h2)+y(j1+1)*h2-y(j1+2)*h)/yy(1) - goto 8 - 7 b0=0.0d0 - 8 b1=b0 - - if (j2 > 1000) stop 'error in subroutine deriv for j2' - - do i=j1,j2 - h=r(i+1)-r(i) - yy(1)=y(i+1)-y(i) - h2=h*h - ha=h-a0 - h2a=h-2.0d0*a0 - h3a=2.0d0*h-3.0d0*a0 - h2b=h2*b0 - s1(i)=h2/ha - s2(i)=-ha/(h2a*h2) - s3(i)=-h*h2a/h3a - f(1,i)=(yy(1)-h*b0)/(h*ha) - f(2,i)=(h2b-yy(1)*(2.0d0*h-a0))/(h*h2*h2a) - f(3,i)=-(h2b-3.0d0*yy(1)*ha)/(h*h3a) - a0=s3(i) - b0=f(3,i) - enddo + ndp = ndis+1 + do 3 nd = 1,ndp + if (nd == 1) then + j1 = 1 + j2 = kdis(1)-2 + else if (nd == ndp) then + j1 = kdis(ndis)+1 + j2 = n-2 + else + j1 = kdis(nd-1)+1 + j2 = kdis(nd)-2 + endif - i=j2+1 - h=r(i+1)-r(i) - yy(1)=y(i+1)-y(i) - h2=h*h - ha=h-a0 - h2a=h*ha - h2b=h2*b0-yy(1)*(2.d0*h-a0) - s1(i)=h2/ha - f(1,i)=(yy(1)-h*b0)/h2a - ha=r(j2)-r(i+1) - yy(1)=-h*ha*(ha+h) - ha=ha*ha - yy(1)=(y(i+1)*(h2-ha)+y(i)*ha-y(j2)*h2)/yy(1) - s3(i)=(yy(1)*h2a+h2b)/(h*h2*(h-2.0d0*a0)) - s13=s1(i)*s3(i) - s2(i)=f(1,i)-s13 - - do j=j1,j2 - k=i-1 - s32=s3(k)*s2(i) - s1(i)=f(3,k)-s32 - s21=s2(k)*s1(i) - s3(k)=f(2,k)-s21 - s13=s1(k)*s3(k) - s2(k)=f(1,k)-s13 - i=k - enddo + if ((j2+1-j1) > 0) goto 11 + + j2 = j2+2 + yy(1) = (y(j2)-y(j1))/(r(j2)-r(j1)) + s1(j1) = yy(1) + s1(j2) = yy(1) + s2(j1) = yy(2) + s2(j2) = yy(2) + s3(j1) = yy(3) + s3(j2) = yy(3) + goto 3 + +11 a0 = 0.0d0 + + if (j1 == 1) then + b0 = 0.0d0 + else + h = r(j1+1)-r(j1) + h2 = r(j1+2)-r(j1) + yy(1) = h*h2*(h2-h) + h = h*h + h2 = h2*h2 + b0 = (y(j1)*(h-h2)+y(j1+1)*h2-y(j1+2)*h)/yy(1) + endif + b1 = b0 + + if (j2 > 1000) stop 'error in subroutine deriv for j2' + + do i = j1,j2 + h = r(i+1)-r(i) + yy(1) = y(i+1)-y(i) + h2 = h*h + ha = h-a0 + h2a = h-2.0d0*a0 + h3a = 2.0d0*h-3.0d0*a0 + h2b = h2*b0 + s1(i) = h2/ha + s2(i) = -ha/(h2a*h2) + s3(i) = -h*h2a/h3a + f(1,i) = (yy(1)-h*b0)/(h*ha) + f(2,i) = (h2b-yy(1)*(2.0d0*h-a0))/(h*h2*h2a) + f(3,i) = -(h2b-3.0d0*yy(1)*ha)/(h*h3a) + a0 = s3(i) + b0 = f(3,i) + enddo + + i = j2+1 + h = r(i+1)-r(i) + yy(1) = y(i+1)-y(i) + h2 = h*h + ha = h-a0 + h2a = h*ha + h2b = h2*b0-yy(1)*(2.d0*h-a0) + s1(i) = h2/ha + f(1,i) = (yy(1)-h*b0)/h2a + ha = r(j2)-r(i+1) + yy(1) = -h*ha*(ha+h) + ha = ha*ha + yy(1) = (y(i+1)*(h2-ha)+y(i)*ha-y(j2)*h2)/yy(1) + s3(i) = (yy(1)*h2a+h2b)/(h*h2*(h-2.0d0*a0)) + s13 = s1(i)*s3(i) + s2(i) = f(1,i)-s13 + + do j = j1,j2 + k = i-1 + s32 = s3(k)*s2(i) + s1(i) = f(3,k)-s32 + s21 = s2(k)*s1(i) + s3(k) = f(2,k)-s21 + s13 = s1(k)*s3(k) + s2(k) = f(1,k)-s13 + i = k + enddo + + s1(i) = b1 + j2 = j2+2 + s1(j2) = yy(1) + s2(j2) = yy(2) + s3(j2) = yy(3) - s1(i)=b1 - j2=j2+2 - s1(j2)=yy(1) - s2(j2)=yy(2) - s3(j2)=yy(3) 3 continue - do i=1,n + do i = 1,n yprime(i)=s1(i) enddo diff --git a/src/specfem3D/pml_compute_accel_contribution.f90 b/src/specfem3D/pml_compute_accel_contribution.f90 index c4627342f..5e140f495 100644 --- a/src/specfem3D/pml_compute_accel_contribution.f90 +++ b/src/specfem3D/pml_compute_accel_contribution.f90 @@ -36,12 +36,13 @@ subroutine pml_compute_accel_contribution_elastic(ispec,ispec_CPML,displ,veloc,r ! Anisotropic-medium PML for vector FETD with modified basis functions, ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006) + use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ + use specfem_par, only: NGLOB_AB,deltat,wgll_cube,jacobian,ibool,rhostore,irregular_element_number,jacobian_regular + use pml_par, only: CPML_regions,d_store_x,d_store_y,d_store_z,K_store_x,K_store_y,K_store_z, & alpha_store_x, alpha_store_y, alpha_store_z, & NSPEC_CPML,accel_elastic_CPML,PML_displ_old,PML_displ_new - use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ,CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, & - CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ implicit none @@ -156,13 +157,14 @@ subroutine pml_compute_accel_contribution_acoustic(ispec,ispec_CPML,potential_ac ! Anisotropic-medium PML for vector FETD with modified basis functions, ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006) + use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ + use specfem_par, only: NGLOB_AB,deltat,wgll_cube,jacobian,ibool,kappastore,irregular_element_number,jacobian_regular + use pml_par, only: CPML_regions,NSPEC_CPML,d_store_x,d_store_y,d_store_z,K_store_x,K_store_y,K_store_z, & alpha_store_x, alpha_store_y, alpha_store_z, & NSPEC_CPML,potential_dot_dot_acoustic_CPML, & PML_potential_acoustic_old,PML_potential_acoustic_new - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, & - CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ implicit none @@ -170,7 +172,6 @@ subroutine pml_compute_accel_contribution_acoustic(ispec,ispec_CPML,potential_ac real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: potential_acoustic,potential_dot_acoustic real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3) :: rmemory_potential_acoustic - ! local parameters integer :: i,j,k,iglob,CPML_region_local,ispec_irreg @@ -782,15 +783,15 @@ subroutine compute_convolution_coef(bb,deltat,coef0,coef1,coef2) ! permanent factors (avoids divisions which are computationally expensive) real(kind=CUSTOM_REAL),parameter :: ONE_OVER_8 = 0.125_CUSTOM_REAL - real(kind=CUSTOM_REAL),parameter :: ONE_OVER_12 = 1._CUSTOM_REAL / 12._CUSTOM_REAL - real(kind=CUSTOM_REAL),parameter :: ONE_OVER_24 = 1._CUSTOM_REAL / 24._CUSTOM_REAL + !real(kind=CUSTOM_REAL),parameter :: ONE_OVER_12 = 1._CUSTOM_REAL / 12._CUSTOM_REAL + !real(kind=CUSTOM_REAL),parameter :: ONE_OVER_24 = 1._CUSTOM_REAL / 24._CUSTOM_REAL real(kind=CUSTOM_REAL),parameter :: ONE_OVER_48 = 1._CUSTOM_REAL / 48._CUSTOM_REAL real(kind=CUSTOM_REAL),parameter :: ONE_OVER_128 = 0.0078125_CUSTOM_REAL real(kind=CUSTOM_REAL),parameter :: ONE_OVER_384 = 1._CUSTOM_REAL / 384._CUSTOM_REAL - real(kind=CUSTOM_REAL),parameter :: ONE_OVER_960 = 1._CUSTOM_REAL / 960._CUSTOM_REAL - real(kind=CUSTOM_REAL),parameter :: ONE_OVER_1920 = 1._CUSTOM_REAL / 1920._CUSTOM_REAL - real(kind=CUSTOM_REAL),parameter :: SEVEN_OVER_3840 = 7._CUSTOM_REAL / 3840._CUSTOM_REAL - real(kind=CUSTOM_REAL),parameter :: FIVE_OVER_11520 = 5._CUSTOM_REAL/11520._CUSTOM_REAL + !real(kind=CUSTOM_REAL),parameter :: ONE_OVER_960 = 1._CUSTOM_REAL / 960._CUSTOM_REAL + !real(kind=CUSTOM_REAL),parameter :: ONE_OVER_1920 = 1._CUSTOM_REAL / 1920._CUSTOM_REAL + !real(kind=CUSTOM_REAL),parameter :: SEVEN_OVER_3840 = 7._CUSTOM_REAL / 3840._CUSTOM_REAL + !real(kind=CUSTOM_REAL),parameter :: FIVE_OVER_11520 = 5._CUSTOM_REAL/11520._CUSTOM_REAL ! helper variables bbpow2 = bb**2 diff --git a/src/specfem3D/pml_compute_memory_variables.f90 b/src/specfem3D/pml_compute_memory_variables.f90 index ade5e27e2..0b1121825 100644 --- a/src/specfem3D/pml_compute_memory_variables.f90 +++ b/src/specfem3D/pml_compute_memory_variables.f90 @@ -42,6 +42,8 @@ subroutine pml_compute_memory_variables_elastic(ispec,ispec_CPML,tempx1,tempy1,t ! Anisotropic-medium PML for vector FETD with modified basis functions, ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006) + use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,FOUR_THIRDS + use specfem_par, only: deltat,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, & kappastore,mustore,irregular_element_number, & jacobian_regular,xix_regular @@ -57,8 +59,6 @@ subroutine pml_compute_memory_variables_elastic(ispec,ispec_CPML,tempx1,tempy1,t PML_duy_dxl_new, PML_duy_dyl_new, PML_duy_dzl_new, & PML_duz_dxl_new, PML_duz_dyl_new, PML_duz_dzl_new - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,FOUR_THIRDS, & - CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ implicit none @@ -387,6 +387,8 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te ! Anisotropic-medium PML for vector FETD with modified basis functions, ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006) + use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ + use specfem_par, only: xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, & deltat,rhostore,irregular_element_number, & jacobian_regular,xix_regular @@ -395,9 +397,6 @@ subroutine pml_compute_memory_variables_acoustic(ispec,ispec_CPML,temp1,temp2,te d_store_x,d_store_y,d_store_z, & alpha_store_x,alpha_store_y,alpha_store_z - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,FOUR_THIRDS, & - CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ - implicit none integer, intent(in) :: ispec,ispec_CPML @@ -559,11 +558,12 @@ subroutine pml_compute_memory_variables_acoustic_elastic(ispec_CPML,iface,iglob, ! Anisotropic-medium PML for vector FETD with modified basis functions, ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006) + use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ + use specfem_par, only: NGLOB_AB,deltat + use pml_par, only: CPML_regions,k_store_x,k_store_y,k_store_z,d_store_x,d_store_y,d_store_z, & alpha_store_x,alpha_store_y,alpha_store_z,PML_displ_old,PML_displ_new - use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ, & - CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ implicit none @@ -655,12 +655,13 @@ subroutine pml_compute_memory_variables_elastic_acoustic(ispec_CPML,iface,iglob, ! Anisotropic-medium PML for vector FETD with modified basis functions, ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006) + use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ + use specfem_par, only: NGLOB_AB,deltat + use pml_par, only: CPML_regions,k_store_x,k_store_y,k_store_z,d_store_x,d_store_y,d_store_z, & alpha_store_x,alpha_store_y,alpha_store_z, & PML_potential_acoustic_old,PML_potential_acoustic_new - use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ, & - CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ implicit none @@ -771,7 +772,7 @@ subroutine lijk_parameter_computation(deltat,kappa_x,d_x,alpha_x,kappa_y,d_y,alp integer :: CPML_X_ONLY_TEMP,CPML_Y_ONLY_TEMP,CPML_Z_ONLY_TEMP, & CPML_XY_ONLY_TEMP,CPML_XZ_ONLY_TEMP,CPML_YZ_ONLY_TEMP,CPML_XYZ_TEMP - logical,parameter :: FIRST_ORDER_CONVOLUTION = .false. + !logical,parameter :: FIRST_ORDER_CONVOLUTION = .false. if (index_ijk == 123) then CPML_X_ONLY_TEMP = CPML_X_ONLY @@ -1195,7 +1196,7 @@ subroutine lxy_interface_parameter_computation(deltat,kappa_x,d_x,alpha_x,kappa_ integer :: CPML_X_ONLY_TEMP,CPML_Y_ONLY_TEMP,CPML_Z_ONLY_TEMP, & CPML_XY_ONLY_TEMP,CPML_XZ_ONLY_TEMP,CPML_YZ_ONLY_TEMP,CPML_XYZ_TEMP - logical,parameter :: FIRST_ORDER_CONVOLUTION = .false. + !logical,parameter :: FIRST_ORDER_CONVOLUTION = .false. if (index_ijk == 12) then CPML_X_ONLY_TEMP = CPML_X_ONLY diff --git a/src/specfem3D/prepare_gpu.f90 b/src/specfem3D/prepare_gpu.f90 index 8ff9f5af4..4799a9e54 100644 --- a/src/specfem3D/prepare_gpu.f90 +++ b/src/specfem3D/prepare_gpu.f90 @@ -54,7 +54,7 @@ subroutine prepare_GPU() call synchronize_all() ! evaluates memory required - call memory_eval_gpu() + call prepare_memory_eval_gpu() ! prepares general fields on GPU ! add GPU support for the C-PML routines @@ -255,174 +255,182 @@ subroutine prepare_GPU() endif call synchronize_all() - contains + end subroutine prepare_GPU + + +! +!---------------------------------------------------------------------------------------------- +! - subroutine memory_eval_gpu() + subroutine prepare_memory_eval_gpu() - implicit none + use specfem_par + use specfem_par_acoustic + use specfem_par_elastic + use specfem_par_poroelastic - ! local parameters - double precision :: memory_size,memory_size_glob - integer,parameter :: NGLL2 = 25 - integer,parameter :: NGLL3 = 125 - integer,parameter :: NGLL3_PADDED = 128 + implicit none - memory_size = 0.d0 + ! local parameters + double precision :: memory_size,memory_size_glob + integer,parameter :: NGLL2 = 25 + integer,parameter :: NGLL3 = 125 + integer,parameter :: NGLL3_PADDED = 128 + + memory_size = 0.d0 + + ! add size of each set of arrays multiplied by the number of such arrays + ! d_hprime_xx,d_hprimewgll_xx + memory_size = memory_size + 2.d0 * NGLL2 * dble(CUSTOM_REAL) + ! padded xix,..gammaz + memory_size = memory_size + 9.d0 * NGLL3_PADDED * NSPEC_IRREGULAR * dble(CUSTOM_REAL) + ! ibool + irregular_element_number + memory_size = memory_size + (NGLL3+1) * NSPEC_AB * dble(SIZE_INTEGER) + ! d_ibool_interfaces_ext_mesh + memory_size = memory_size + num_interfaces_ext_mesh * max_nibool_interfaces_ext_mesh * dble(SIZE_INTEGER) + ! ispec_is_inner + memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER) + + if (STACEY_ABSORBING_CONDITIONS) then + ! d_abs_boundary_ispec + memory_size = memory_size + num_abs_boundary_faces * dble(SIZE_INTEGER) + ! d_abs_boundary_ijk + memory_size = memory_size + num_abs_boundary_faces * 3.d0 * NGLL2 * dble(SIZE_INTEGER) + ! d_abs_boundary_normal + memory_size = memory_size + num_abs_boundary_faces * NDIM * NGLL2 * dble(CUSTOM_REAL) + ! d_abs_boundary_jacobian2Dw + memory_size = memory_size + num_abs_boundary_faces * NGLL2 * dble(CUSTOM_REAL) + endif - ! add size of each set of arrays multiplied by the number of such arrays - ! d_hprime_xx,d_hprimewgll_xx - memory_size = memory_size + 2.d0 * NGLL2 * dble(CUSTOM_REAL) - ! padded xix,..gammaz - memory_size = memory_size + 9.d0 * NGLL3_PADDED * NSPEC_IRREGULAR * dble(CUSTOM_REAL) - ! ibool + irregular_element_number - memory_size = memory_size + (NGLL3+1) * NSPEC_AB * dble(SIZE_INTEGER) - ! d_ibool_interfaces_ext_mesh - memory_size = memory_size + num_interfaces_ext_mesh * max_nibool_interfaces_ext_mesh * dble(SIZE_INTEGER) - ! ispec_is_inner + ! sources + ! d_sourcearrays + memory_size = memory_size + NGLL3 * NSOURCES * NDIM * dble(CUSTOM_REAL) + ! d_islice_selected_source,d_ispec_selected_source + memory_size = memory_size + 2.0 * NSOURCES * dble(SIZE_INTEGER) + + ! receivers + ! d_ispec_selected_rec_loc + memory_size = memory_size + nrec_local * dble(SIZE_INTEGER) + ! d_ispec_selected_rec + memory_size = memory_size + nrec * dble(SIZE_INTEGER) + + ! d_hxir, d_hetar, d_hgammar + memory_size = memory_size + 3.d0 * NGLLX * nrec_local * dble(CUSTOM_REAL) + ! d_nu + memory_size = memory_size + NDIM * NDIM * nrec_local * dble(CUSTOM_REAL) + ! d_ispec_selected_rec_loc + memory_size = memory_size + nrec_local * dble(SIZE_INTEGER) + ! d_seismograms_d,d_seismograms_v,d_seismograms_a,d_seismograms_p + if (SAVE_SEISMOGRAMS_DISPLACEMENT) & + memory_size = memory_size + NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) + if (SAVE_SEISMOGRAMS_VELOCITY) & + memory_size = memory_size + NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) + if (SAVE_SEISMOGRAMS_ACCELERATION) & + memory_size = memory_size + NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) + if (SAVE_SEISMOGRAMS_PRESSURE) & + memory_size = memory_size + NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) * NB_RUNS_ACOUSTIC_GPU + + ! acoustic simulations + if (ACOUSTIC_SIMULATION) then + ! d_potential_acoustic,d_potential_dot_acoustic,d_potential_dot_dot_acoustic + memory_size = memory_size + 3.d0 * NGLOB_AB * dble(CUSTOM_REAL) * NB_RUNS_ACOUSTIC_GPU + ! d_rmass_acoustic + memory_size = memory_size + NGLOB_AB * dble(CUSTOM_REAL) + ! padded d_rhostore + memory_size = memory_size + NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL) + ! d_kappastore + memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) + ! d_phase_ispec_inner_acoustic + memory_size = memory_size + 2.d0 * num_phase_ispec_acoustic * dble(SIZE_INTEGER) + ! d_ispec_is_acoustic memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER) - if (STACEY_ABSORBING_CONDITIONS) then - ! d_abs_boundary_ispec - memory_size = memory_size + num_abs_boundary_faces * dble(SIZE_INTEGER) - ! d_abs_boundary_ijk - memory_size = memory_size + num_abs_boundary_faces * 3.d0 * NGLL2 * dble(SIZE_INTEGER) - ! d_abs_boundary_normal - memory_size = memory_size + num_abs_boundary_faces * NDIM * NGLL2 * dble(CUSTOM_REAL) - ! d_abs_boundary_jacobian2Dw - memory_size = memory_size + num_abs_boundary_faces * NGLL2 * dble(CUSTOM_REAL) - endif + endif - ! sources - ! d_sourcearrays - memory_size = memory_size + NGLL3 * NSOURCES * NDIM * dble(CUSTOM_REAL) - ! d_islice_selected_source,d_ispec_selected_source - memory_size = memory_size + 2.0 * NSOURCES * dble(SIZE_INTEGER) - - ! receivers - ! d_ispec_selected_rec_loc - memory_size = memory_size + nrec_local * dble(SIZE_INTEGER) - ! d_ispec_selected_rec - memory_size = memory_size + nrec * dble(SIZE_INTEGER) - - ! d_hxir, d_hetar, d_hgammar - memory_size = memory_size + 3.d0 * NGLLX * nrec_local * dble(CUSTOM_REAL) - ! d_nu - memory_size = memory_size + NDIM * NDIM * nrec_local * dble(CUSTOM_REAL) - ! d_ispec_selected_rec_loc - memory_size = memory_size + nrec_local * dble(SIZE_INTEGER) - ! d_seismograms_d,d_seismograms_v,d_seismograms_a,d_seismograms_p - if (SAVE_SEISMOGRAMS_DISPLACEMENT) & - memory_size = memory_size + NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) - if (SAVE_SEISMOGRAMS_VELOCITY) & - memory_size = memory_size + NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) - if (SAVE_SEISMOGRAMS_ACCELERATION) & - memory_size = memory_size + NDIM * NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) - if (SAVE_SEISMOGRAMS_PRESSURE) & - memory_size = memory_size + NTSTEP_BETWEEN_OUTPUT_SEISMOS * nrec_local * dble(CUSTOM_REAL) * NB_RUNS_ACOUSTIC_GPU - - ! acoustic simulations - if (ACOUSTIC_SIMULATION) then - ! d_potential_acoustic,d_potential_dot_acoustic,d_potential_dot_dot_acoustic - memory_size = memory_size + 3.d0 * NGLOB_AB * dble(CUSTOM_REAL) * NB_RUNS_ACOUSTIC_GPU - ! d_rmass_acoustic - memory_size = memory_size + NGLOB_AB * dble(CUSTOM_REAL) - ! padded d_rhostore - memory_size = memory_size + NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL) - ! d_kappastore - memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) - ! d_phase_ispec_inner_acoustic - memory_size = memory_size + 2.d0 * num_phase_ispec_acoustic * dble(SIZE_INTEGER) - ! d_ispec_is_acoustic - memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER) + ! elastic simulations + if (ELASTIC_SIMULATION) then + ! d_displ,.. + memory_size = memory_size + 3.d0 * NDIM * NGLOB_AB * dble(CUSTOM_REAL) + ! d_send_accel_buffer + memory_size = memory_size + 3.d0 * num_interfaces_ext_mesh * max_nibool_interfaces_ext_mesh * dble(CUSTOM_REAL) + ! d_rmassx,.. + memory_size = memory_size + 3.d0 * NGLOB_AB * dble(CUSTOM_REAL) + ! d_ispec_is_elastic + memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER) + ! d_phase_ispec_inner_elastic + memory_size = memory_size + 2.d0 * num_phase_ispec_elastic * dble(SIZE_INTEGER) + if (STACEY_ABSORBING_CONDITIONS .or. PML_CONDITIONS) then + ! d_rho_vp,.. + memory_size = memory_size + 2.d0 * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) endif - ! elastic simulations - if (ELASTIC_SIMULATION) then - ! d_displ,.. - memory_size = memory_size + 3.d0 * NDIM * NGLOB_AB * dble(CUSTOM_REAL) - ! d_send_accel_buffer - memory_size = memory_size + 3.d0 * num_interfaces_ext_mesh * max_nibool_interfaces_ext_mesh * dble(CUSTOM_REAL) - ! d_rmassx,.. - memory_size = memory_size + 3.d0 * NGLOB_AB * dble(CUSTOM_REAL) - ! d_ispec_is_elastic - memory_size = memory_size + NSPEC_AB * dble(SIZE_INTEGER) - ! d_phase_ispec_inner_elastic - memory_size = memory_size + 2.d0 * num_phase_ispec_elastic * dble(SIZE_INTEGER) - - if (STACEY_ABSORBING_CONDITIONS .or. PML_CONDITIONS) then - ! d_rho_vp,.. - memory_size = memory_size + 2.d0 * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) - endif - - ! padded kappav,muv - memory_size = memory_size + 2.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL) - - if (COMPUTE_AND_STORE_STRAIN) then - ! d_epsilondev_xx,.. - memory_size = memory_size + 5.d0 * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) - endif - if (ATTENUATION) then - ! d_R_xx,.. - memory_size = memory_size + 5.d0 * size(R_xx) * dble(CUSTOM_REAL) - ! d_factor_common - memory_size = memory_size + N_SLS * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) - ! alphaval,.. - memory_size = memory_size + 3.d0 * N_SLS * dble(CUSTOM_REAL) - endif - if (ANISOTROPY) then - ! padded d_c11store,.. - memory_size = memory_size + 21.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL) - endif - if (APPROXIMATE_OCEAN_LOAD) then - ! d_rmass_ocean_load - memory_size = memory_size + NGLOB_AB * dble(CUSTOM_REAL) - ! d_free_surface_normal - memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL) - endif - endif + ! padded kappav,muv + memory_size = memory_size + 2.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL) - ! noise simulations - if (NOISE_TOMOGRAPHY > 0) then - ! d_free_surface_ispec - memory_size = memory_size + num_free_surface_faces * dble(SIZE_INTEGER) - ! d_free_surface_ijk - memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(SIZE_INTEGER) - ! d_noise_surface_movie + if (COMPUTE_AND_STORE_STRAIN) then + ! d_epsilondev_xx,.. + memory_size = memory_size + 5.d0 * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) + endif + if (ATTENUATION) then + ! d_R_xx,.. + memory_size = memory_size + 5.d0 * size(R_xx) * dble(CUSTOM_REAL) + ! d_factor_common + memory_size = memory_size + N_SLS * NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) + ! alphaval,.. + memory_size = memory_size + 3.d0 * N_SLS * dble(CUSTOM_REAL) + endif + if (ANISOTROPY) then + ! padded d_c11store,.. + memory_size = memory_size + 21.d0 * NGLL3_PADDED * NSPEC_AB * dble(CUSTOM_REAL) + endif + if (APPROXIMATE_OCEAN_LOAD) then + ! d_rmass_ocean_load + memory_size = memory_size + NGLOB_AB * dble(CUSTOM_REAL) + ! d_free_surface_normal memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL) - if (NOISE_TOMOGRAPHY == 1) then - ! d_noise_sourcearray - memory_size = memory_size + 3.d0 * NGLL3 * NSTEP * dble(CUSTOM_REAL) - endif - if (NOISE_TOMOGRAPHY > 1) then - ! d_normal_x_noise,.. - memory_size = memory_size + 5.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL) - endif - if (NOISE_TOMOGRAPHY == 3) then - ! d_sigma_kl - memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) - endif endif + endif - if (GRAVITY) then - ! d_minus_deriv_gravity,d_minus_g - memory_size = memory_size + 2.d0 * NGLOB_AB * dble(CUSTOM_REAL) + ! noise simulations + if (NOISE_TOMOGRAPHY > 0) then + ! d_free_surface_ispec + memory_size = memory_size + num_free_surface_faces * dble(SIZE_INTEGER) + ! d_free_surface_ijk + memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(SIZE_INTEGER) + ! d_noise_surface_movie + memory_size = memory_size + 3.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL) + if (NOISE_TOMOGRAPHY == 1) then + ! d_noise_sourcearray + memory_size = memory_size + 3.d0 * NGLL3 * NSTEP * dble(CUSTOM_REAL) + endif + if (NOISE_TOMOGRAPHY > 1) then + ! d_normal_x_noise,.. + memory_size = memory_size + 5.d0 * NGLL2 * num_free_surface_faces * dble(CUSTOM_REAL) endif + if (NOISE_TOMOGRAPHY == 3) then + ! d_sigma_kl + memory_size = memory_size + NGLL3 * NSPEC_AB * dble(CUSTOM_REAL) + endif + endif - ! poor estimate for kernel simulations... - if (SIMULATION_TYPE == 3) memory_size = 2.d0 * memory_size + if (GRAVITY) then + ! d_minus_deriv_gravity,d_minus_g + memory_size = memory_size + 2.d0 * NGLOB_AB * dble(CUSTOM_REAL) + endif - ! maximum of all processes (may vary e.g. due to different nrec_local, num_abs_boundary_faces, ..) - call max_all_dp(memory_size,memory_size_glob) + ! poor estimate for kernel simulations... + if (SIMULATION_TYPE == 3) memory_size = 2.d0 * memory_size - ! user output - if (myrank == 0) then - write(IMAIN,*) - write(IMAIN,*) ' minimum memory requested : ',memory_size_glob / 1024. / 1024.,'MB per process' - write(IMAIN,*) - call flush_IMAIN() - endif + ! maximum of all processes (may vary e.g. due to different nrec_local, num_abs_boundary_faces, ..) + call max_all_dp(memory_size,memory_size_glob) - end subroutine memory_eval_gpu + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) ' minimum memory requested : ',memory_size_glob / 1024. / 1024.,'MB per process' + write(IMAIN,*) + call flush_IMAIN() + endif - end subroutine prepare_GPU + end subroutine prepare_memory_eval_gpu diff --git a/src/specfem3D/prepare_timerun.F90 b/src/specfem3D/prepare_timerun.F90 index c01e4498b..c32f2531d 100644 --- a/src/specfem3D/prepare_timerun.F90 +++ b/src/specfem3D/prepare_timerun.F90 @@ -556,7 +556,7 @@ end subroutine prepare_timerun_lddrk subroutine prepare_timerun_pml() - use constants, only: IMAIN,NGNOD_EIGHT_CORNERS + use constants, only: IMAIN use specfem_par, only: myrank,SIMULATION_TYPE,GPU_MODE,UNDO_ATTENUATION_AND_OR_PML,PML_CONDITIONS use pml_par diff --git a/src/specfem3D/read_stations.f90 b/src/specfem3D/read_stations.f90 index 822184a80..1f2baa2f5 100644 --- a/src/specfem3D/read_stations.f90 +++ b/src/specfem3D/read_stations.f90 @@ -132,7 +132,7 @@ subroutine read_stations_SU_from_previous_run(nrec,station_name,network_name, & xi_receiver,eta_receiver,gamma_receiver, & nu,is_done_stations) - use constants, only: NDIM,IIN,IOUT_SU, & + use constants, only: NDIM,IOUT_SU, & MAX_LENGTH_STATION_NAME,MAX_LENGTH_NETWORK_NAME,IMAIN,OUTPUT_FILES,MAX_STRING_LEN use specfem_par, only: myrank diff --git a/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 b/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 index 5d236e65a..3a62e64e6 100644 --- a/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 +++ b/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 @@ -59,8 +59,8 @@ program smooth_sem use constants, only: USE_QUADRATURE_RULE_FOR_SMOOTHING - use postprocess_par, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,NGLLSQUARE, & - MAX_STRING_LEN,IIN,IOUT,GAUSSALPHA,GAUSSBETA,PI,TWO_PI,MAX_KERNEL_NAMES + use postprocess_par, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ, & + MAX_STRING_LEN,IIN,IOUT,GAUSSALPHA,GAUSSBETA,PI,MAX_KERNEL_NAMES use specfem_par use specfem_par_elastic, only: ELASTIC_SIMULATION,ispec_is_elastic,rho_vp,rho_vs,min_resolved_period diff --git a/src/tomography/read_parameters_tomo.f90 b/src/tomography/read_parameters_tomo.f90 index e323284c8..967fe9ab5 100644 --- a/src/tomography/read_parameters_tomo.f90 +++ b/src/tomography/read_parameters_tomo.f90 @@ -42,13 +42,13 @@ subroutine read_parameters_tomo() call get_command_argument(1,s_step_fac) if (trim(s_step_fac) == '') then - call usage() + call tomo_usage() endif ! read in parameter information read(s_step_fac,*,iostat=ier) step_fac if (ier /= 0) then - call usage() + call tomo_usage() endif ! safety check @@ -96,9 +96,15 @@ subroutine read_parameters_tomo() close(IOUT) endif -contains + end subroutine read_parameters_tomo - subroutine usage() +! +!------------------------------------------------------------------------------------------------- +! + + subroutine tomo_usage() + + use tomography_par implicit none @@ -116,7 +122,6 @@ subroutine usage() call synchronize_all() call exit_MPI(myrank,'Error usage: add_model step_factor') - end subroutine usage + end subroutine tomo_usage -end subroutine read_parameters_tomo