[cig-commits] r17959 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh_SCOTCH generate_databases shared specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Feb 23 13:22:43 PST 2011
Author: danielpeter
Date: 2011-02-23 13:22:42 -0800 (Wed, 23 Feb 2011)
New Revision: 17959
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90
seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90
seismo/3D/SPECFEM3D/trunk/src/shared/sort_array_coordinates.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
Log:
updates minimum period calculation for SPECFEM3D/
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -564,12 +564,12 @@
nb_edges = xadj(nspec+1)
- ! allocates & initializes partioning of elements
+ ! allocates & initializes partioning of elements
allocate(part(1:nspec))
part(:) = -1
- ! initializes
- ! elements load array
+ ! initializes
+ ! elements load array
allocate(elmnts_load(1:nspec))
! uniform load
@@ -579,7 +579,7 @@
call acoustic_elastic_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
mat(1,:),mat_prop,undef_mat_prop)
- ! SCOTCH partitioning
+ ! SCOTCH partitioning
! we use default strategy for partitioning, thus omit specifing explicit strategy .
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -90,20 +90,23 @@
if ( .not.is_neighbour ) then
if ( adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour+m) == nodes_elmnts(l+j*nsize) ) then
is_neighbour = .true.
-
end if
end if
end do
if ( .not.is_neighbour ) then
- adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
+ adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour &
+ + xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
- xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
- if (xadj(nodes_elmnts(k+j*nsize))>sup_neighbour) stop 'ERROR : too much neighbours per element, modify the mesh.'
+ xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
+ if (xadj(nodes_elmnts(k+j*nsize))>sup_neighbour) &
+ stop 'ERROR : too much neighbours per element, modify the mesh.'
- adjncy(nodes_elmnts(l+j*nsize)*sup_neighbour+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
+ adjncy(nodes_elmnts(l+j*nsize)*sup_neighbour &
+ + xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
- if (xadj(nodes_elmnts(l+j*nsize))>sup_neighbour) stop 'ERROR : too much neighbours per element, modify the mesh.'
+ if (xadj(nodes_elmnts(l+j*nsize))>sup_neighbour) &
+ stop 'ERROR : too much neighbours per element, modify the mesh.'
end if
end if
end do
@@ -144,7 +147,7 @@
! allocates local numbering array
allocate(glob2loc_elmnts(0:nelmnts-1))
- ! initializes number of local points per partition
+ ! initializes number of local elements per partition
do num_part = 0, nparts-1
num_loc(num_part) = 0
end do
@@ -196,14 +199,12 @@
glob2loc_nodes_nparts(num_node) = size_glob2loc_nodes
do el = 0, nnodes_elmnts(num_node)-1
parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
-
end do
do num_part = 0, nparts-1
if ( parts_node(num_part) == 1 ) then
size_glob2loc_nodes = size_glob2loc_nodes + 1
parts_node(num_part) = 0
-
end if
end do
@@ -224,7 +225,6 @@
do num_node = 0, nnodes-1
do el = 0, nnodes_elmnts(num_node)-1
parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
-
end do
do num_part = 0, nparts-1
@@ -332,7 +332,8 @@
do num_node = 0, esize-1
do num_node_bis = 0, esize-1
if ( elmnts(el*esize+num_node) == elmnts(adjncy(el_adj)*esize+num_node_bis) ) then
- tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+3+ncommon_nodes) &
+ tab_interfaces(tab_size_interfaces(num_interface)*7 &
+ +num_edge*7+3+ncommon_nodes) &
= elmnts(el*esize+num_node)
ncommon_nodes = ncommon_nodes + 1
end if
@@ -432,8 +433,10 @@
else
is_acoustic_el_adj = .false.
end if
- ! adds element if neighbor element has same material acoustic/not-acoustic and lies in next partition
- if ( (part(adjncy(el_adj)) == num_part_bis) .and. (is_acoustic_el .eqv. is_acoustic_el_adj) ) then
+ ! adds element if neighbor element has same material acoustic/not-acoustic
+ ! and lies in next partition
+ if ( (part(adjncy(el_adj)) == num_part_bis) .and. &
+ (is_acoustic_el .eqv. is_acoustic_el_adj) ) then
num_edge = num_edge + 1
end if
end do
@@ -478,14 +481,16 @@
else
is_acoustic_el_adj = .false.
end if
- if ( (part(adjncy(el_adj)) == num_part_bis) .and. (is_acoustic_el .eqv. is_acoustic_el_adj) ) then
+ if ( (part(adjncy(el_adj)) == num_part_bis) .and. &
+ (is_acoustic_el .eqv. is_acoustic_el_adj) ) then
tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+0) = el
tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+1) = adjncy(el_adj)
ncommon_nodes = 0
do num_node = 0, esize-1
do num_node_bis = 0, esize-1
if ( elmnts(el*esize+num_node) == elmnts(adjncy(el_adj)*esize+num_node_bis) ) then
- tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+3+ncommon_nodes) &
+ tab_interfaces(tab_size_interfaces(num_interface)*7 &
+ +num_edge*7+3+ncommon_nodes) &
= elmnts(el*esize+num_node)
ncommon_nodes = ncommon_nodes + 1
end if
@@ -546,7 +551,8 @@
do i = 0, nnodes-1
do j = glob2loc_nodes_nparts(i), glob2loc_nodes_nparts(i+1)-1
if ( glob2loc_nodes_parts(j) == iproc ) then
- write(IIN_database,*) glob2loc_nodes(j)+1, nodes_coords(1,i+1), nodes_coords(2,i+1), nodes_coords(3,i+1)
+ write(IIN_database,*) glob2loc_nodes(j)+1, nodes_coords(1,i+1), &
+ nodes_coords(2,i+1), nodes_coords(3,i+1)
end if
end do
end do
@@ -558,7 +564,8 @@
!--------------------------------------------------
! Write material properties in the Database
!--------------------------------------------------
- subroutine write_material_props_database(IIN_database,count_def_mat,count_undef_mat, mat_prop, undef_mat_prop)
+ subroutine write_material_props_database(IIN_database,count_def_mat, &
+ count_undef_mat, mat_prop, undef_mat_prop)
integer, intent(in) :: IIN_database
integer, intent(in) :: count_def_mat,count_undef_mat
@@ -587,7 +594,8 @@
!--------------------------------------------------
- ! Write elements on boundaries (and their four nodes on boundaries) pertaining to iproc partition in the corresponding Database
+ ! Write elements on boundaries (and their four nodes on boundaries)
+ ! pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
subroutine write_boundaries_database(IIN_database, iproc, nelmnts, nspec2D_xmin, nspec2D_xmax, &
nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top, &
@@ -601,7 +609,8 @@
integer, intent(in) :: IIN_database
integer, intent(in) :: iproc
integer(long), intent(in) :: nelmnts
- integer, intent(in) :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top
+ integer, intent(in) :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, &
+ nspec2D_ymax, nspec2D_bottom, nspec2D_top
integer, dimension(nspec2D_xmin), intent(in) :: ibelm_xmin
integer, dimension(nspec2D_xmax), intent(in) :: ibelm_xmax
integer, dimension(nspec2D_ymin), intent(in) :: ibelm_ymin
@@ -678,126 +687,150 @@
! we need to have the arg of glob2loc_elmnts start at 0 ==> glob2loc_nodes(ibelm_** -1 )
do i=1,nspec2D_xmin
if(part(ibelm_xmin(i)) == iproc) then
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node1 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node2 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node3 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- write(IIN_database,*) glob2loc_elmnts(ibelm_xmin(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
+ write(IIN_database,*) glob2loc_elmnts(ibelm_xmin(i)-1)+1, &
+ loc_node1, loc_node2, loc_node3, loc_node4
end if
end do
do i=1,nspec2D_xmax
if(part(ibelm_xmax(i)) == iproc) then
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node1 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node2 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node3 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- write(IIN_database,*) glob2loc_elmnts(ibelm_xmax(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
+ write(IIN_database,*) glob2loc_elmnts(ibelm_xmax(i)-1)+1, &
+ loc_node1, loc_node2, loc_node3, loc_node4
end if
end do
do i=1,nspec2D_ymin
if(part(ibelm_ymin(i)) == iproc) then
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node1 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node2 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node3 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- write(IIN_database,*) glob2loc_elmnts(ibelm_ymin(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
+ write(IIN_database,*) glob2loc_elmnts(ibelm_ymin(i)-1)+1, &
+ loc_node1, loc_node2, loc_node3, loc_node4
end if
end do
do i=1,nspec2D_ymax
if(part(ibelm_ymax(i)) == iproc) then
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node1 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node2 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node3 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- write(IIN_database,*) glob2loc_elmnts(ibelm_ymax(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
+ write(IIN_database,*) glob2loc_elmnts(ibelm_ymax(i)-1)+1, &
+ loc_node1, loc_node2, loc_node3, loc_node4
end if
end do
do i=1,nspec2D_bottom
if(part(ibelm_bottom(i)) == iproc) then
- do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node1 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node2 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node3 = glob2loc_nodes(j)+1
end if
end do
- do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i))-1
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i)-1), &
+ glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i))-1
if (glob2loc_nodes_parts(j) == iproc ) then
loc_node4 = glob2loc_nodes(j)+1
end if
@@ -888,7 +921,8 @@
! format:
! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
- write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(1,i+1), num_modele(2,i+1),(loc_nodes(k)+1, k=0,ngnod-1)
+ write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(1,i+1), &
+ num_modele(2,i+1),(loc_nodes(k)+1, k=0,ngnod-1)
end if
end do
end if
@@ -901,9 +935,11 @@
!--------------------------------------------------
! Write interfaces (element and common nodes) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
- subroutine write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, iproc, ninterfaces, &
- my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
- glob2loc_nodes, num_phase, nparts)
+ subroutine write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, &
+ iproc, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, num_phase, nparts)
! include './constants_decompose_mesh_SCOTCH.h'
@@ -945,7 +981,8 @@
! sets flag
my_interfaces(num_interface) = 1
! sets number of elements on interface
- my_nb_interfaces(num_interface) = tab_size_interfaces(num_interface+1) - tab_size_interfaces(num_interface)
+ my_nb_interfaces(num_interface) = tab_size_interfaces(num_interface+1) &
+ - tab_size_interfaces(num_interface)
end if
num_interface = num_interface + 1
end do
@@ -1008,7 +1045,8 @@
local_nodes(1) = glob2loc_nodes(l)+1
end if
end do
- write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), local_nodes(1), -1, -1, -1
+ write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), &
+ local_nodes(1), -1, -1, -1
case (2)
! edge element
do l = glob2loc_nodes_nparts(tab_interfaces(k*7+3)), &
@@ -1023,7 +1061,8 @@
local_nodes(2) = glob2loc_nodes(l)+1
end if
end do
- write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), local_nodes(1), local_nodes(2), -1, -1
+ write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), &
+ local_nodes(1), local_nodes(2), -1, -1
case (4)
! face element
count_faces = count_faces + 1
@@ -1142,7 +1181,8 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- write(IIN_database,*) glob2loc_elmnts(ibelm_moho(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
+ write(IIN_database,*) glob2loc_elmnts(ibelm_moho(i)-1)+1, &
+ loc_node1, loc_node2, loc_node3, loc_node4
end if
end do
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -25,12 +25,12 @@
!=====================================================================
subroutine get_MPI(myrank,nglob,nspec,ibool, &
- nelmnts_ext_mesh,elmnts_ext_mesh, &
- my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
- ibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh, &
- num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
- my_neighbours_ext_mesh,NPROC)
+ nelmnts_ext_mesh,elmnts_ext_mesh, &
+ my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh, &
+ num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,NPROC)
! sets up the MPI interface for communication between partitions
@@ -49,18 +49,17 @@
integer :: num_interfaces_ext_mesh,max_interface_size_ext_mesh
integer, dimension(num_interfaces_ext_mesh) :: my_nelmnts_neighbours_ext_mesh
- integer, dimension(6,max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: my_interfaces_ext_mesh
+ integer, dimension(6,max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: &
+ my_interfaces_ext_mesh
integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
- integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: &
+ ibool_interfaces_ext_mesh
- !integer :: nnodes_ext_mesh
- !double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
-
-!local parameters
+ !local parameters
double precision, dimension(:), allocatable :: xp,yp,zp
double precision, dimension(:), allocatable :: work_ext_mesh
@@ -68,8 +67,8 @@
integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh_true
! for MPI buffers
- integer, dimension(:), allocatable :: reorder_interface_ext_mesh,ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
- integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
+ integer, dimension(:), allocatable :: reorder_interface_ext_mesh, &
+ ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
logical, dimension(:), allocatable :: ifseg
integer :: iinterface,ilocnum
integer :: num_points1, num_points2
@@ -81,8 +80,10 @@
real(kind=CUSTOM_REAL), dimension(:),allocatable :: test_flag_cr
integer, dimension(:,:), allocatable :: ibool_interfaces_dummy
-! gets global indices for points on MPI interfaces (defined by my_interfaces_ext_mesh) between different partitions
-! and stores them in ibool_interfaces_ext_mesh & nibool_interfaces_ext_mesh (number of total points)
+ ! gets global indices for points on MPI interfaces
+ ! (defined by my_interfaces_ext_mesh) between different partitions
+ ! and stores them in ibool_interfaces_ext_mesh & nibool_interfaces_ext_mesh
+ ! (number of total points)
call prepare_assemble_MPI( nelmnts_ext_mesh,elmnts_ext_mesh, &
ibool,nglob,ESIZE, &
num_interfaces_ext_mesh, max_interface_size_ext_mesh, &
@@ -92,7 +93,7 @@
allocate(nibool_interfaces_ext_mesh_true(num_interfaces_ext_mesh))
-! sorts ibool comm buffers lexicographically for all MPI interfaces
+ ! sorts ibool comm buffers lexicographically for all MPI interfaces
num_points1 = 0
num_points2 = 0
do iinterface = 1, num_interfaces_ext_mesh
@@ -103,7 +104,6 @@
allocate(locval(nibool_interfaces_ext_mesh(iinterface)))
allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)))
allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ibool_interface_ext_mesh_dummy(nibool_interfaces_ext_mesh(iinterface)))
allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
@@ -120,7 +120,8 @@
! of total number of points nibool_interfaces_ext_mesh_true(iinterface)
call sort_array_coordinates(nibool_interfaces_ext_mesh(iinterface),xp,yp,zp, &
ibool_interfaces_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
- reorder_interface_ext_mesh,locval,ifseg,nibool_interfaces_ext_mesh_true(iinterface), &
+ reorder_interface_ext_mesh,locval,ifseg, &
+ nibool_interfaces_ext_mesh_true(iinterface), &
ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh)
! checks that number of MPI points are still the same
@@ -140,7 +141,6 @@
deallocate(locval)
deallocate(ifseg)
deallocate(reorder_interface_ext_mesh)
- deallocate(ibool_interface_ext_mesh_dummy)
deallocate(ind_ext_mesh)
deallocate(ninseg_ext_mesh)
deallocate(iwork_ext_mesh)
@@ -157,7 +157,7 @@
write(IMAIN,*) ' total MPI interface points: ',ilocnum
endif
-! checks with assembly of test fields
+ ! checks with assembly of test fields
allocate(test_flag(nglob),test_flag_cr(nglob))
test_flag(:) = 0
test_flag_cr(:) = 0._CUSTOM_REAL
@@ -189,9 +189,11 @@
count = 0
do iinterface = 1, num_interfaces_ext_mesh
- ibool_interfaces_dummy(:,iinterface) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,iinterface)
+ ibool_interfaces_dummy(:,iinterface) = &
+ ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,iinterface)
count = count + nibool_interfaces_ext_mesh(iinterface)
- !write(*,*) myrank,'interfaces ',iinterface,nibool_interfaces_ext_mesh(iinterface),max_nibool_interfaces_ext_mesh
+ !write(*,*) myrank,'interfaces ',iinterface, &
+ ! nibool_interfaces_ext_mesh(iinterface),max_nibool_interfaces_ext_mesh
enddo
call sync_all()
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -85,12 +85,14 @@
! send messages
do iinterface = 1, num_interfaces_ext_mesh
+ ! non-blocking synchronous send request
call issend_cr(buffer_send_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
nibool_interfaces_ext_mesh(iinterface), &
my_neighbours_ext_mesh(iinterface), &
itag, &
request_send_scalar_ext_mesh(iinterface) &
)
+ ! receive request
call irecv_cr(buffer_recv_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
nibool_interfaces_ext_mesh(iinterface), &
my_neighbours_ext_mesh(iinterface), &
@@ -99,7 +101,7 @@
)
enddo
- ! wait for communications completion
+ ! wait for communications completion (recv)
do iinterface = 1, num_interfaces_ext_mesh
call wait_req(request_recv_scalar_ext_mesh(iinterface))
enddo
@@ -169,18 +171,21 @@
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
- buffer_send_scalar_ext_mesh(ipoin,iinterface) = array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
+ buffer_send_scalar_ext_mesh(ipoin,iinterface) = &
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
enddo
enddo
! send messages
do iinterface = 1, num_interfaces_ext_mesh
+ ! non-blocking synchronous send request
call issend_i(buffer_send_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
nibool_interfaces_ext_mesh(iinterface), &
my_neighbours_ext_mesh(iinterface), &
itag, &
request_send_scalar_ext_mesh(iinterface) &
)
+ ! receive request
call irecv_i(buffer_recv_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
nibool_interfaces_ext_mesh(iinterface), &
my_neighbours_ext_mesh(iinterface), &
@@ -198,7 +203,8 @@
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
- array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) &
+ + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
enddo
enddo
@@ -256,18 +262,21 @@
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
- buffer_send_scalar_ext_mesh(ipoin,iinterface) = array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
+ buffer_send_scalar_ext_mesh(ipoin,iinterface) = &
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
enddo
enddo
! send messages
do iinterface = 1, num_interfaces_ext_mesh
+ ! non-blocking synchronous send request
call issend_cr(buffer_send_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
nibool_interfaces_ext_mesh(iinterface), &
my_neighbours_ext_mesh(iinterface), &
itag, &
request_send_scalar_ext_mesh(iinterface) &
)
+ ! receive request
call irecv_cr(buffer_recv_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
nibool_interfaces_ext_mesh(iinterface), &
my_neighbours_ext_mesh(iinterface), &
@@ -316,7 +325,7 @@
! assemble only if more than one partition
if(NPROC > 1) then
- ! wait for communications completion
+ ! wait for communications completion (recv)
do iinterface = 1, num_interfaces_ext_mesh
call wait_req(request_recv_scalar_ext_mesh(iinterface))
enddo
@@ -325,7 +334,8 @@
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
- array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) &
+ + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -45,26 +45,34 @@
real(kind=CUSTOM_REAL) :: model_speed_max,min_resolved_period
! local parameters
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ):: vp_elem,vs_elem
- real(kind=CUSTOM_REAL), dimension(1) :: val_min,val_max
real(kind=CUSTOM_REAL) :: vpmin,vpmax,vsmin,vsmax,vpmin_glob,vpmax_glob,vsmin_glob,vsmax_glob
- real(kind=CUSTOM_REAL) :: distance_min,distance_max,distance_min_glob,distance_max_glob,dx !,dy,dz
+ real(kind=CUSTOM_REAL) :: distance_min,distance_max,distance_min_glob,distance_max_glob
+ real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max,elemsize_min_glob,elemsize_max_glob
real(kind=CUSTOM_REAL) :: cmax,cmax_glob,pmax,pmax_glob
- real(kind=CUSTOM_REAL) :: dt_suggested,dt_suggested_glob
+ real(kind=CUSTOM_REAL) :: dt_suggested,dt_suggested_glob,avg_distance
logical:: DT_PRESENT
integer :: myrank
integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum
integer :: NGLOB_AB_global_min,NGLOB_AB_global_max,NGLOB_AB_global_sum
- integer :: i,j,k,ii,jj,kk,ispec,iglob_a,iglob_b,sizeprocs
+ integer :: ispec,sizeprocs
-! estimation of time step and period resolved
- real(kind=CUSTOM_REAL),parameter :: COURANT_SUGGESTED = 0.3
- real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
+ !********************************************************************************
+
+ ! empirical choice for distorted elements to estimate time step and period resolved:
+ ! courant number for time step estimate
+ real(kind=CUSTOM_REAL),parameter :: COURANT_SUGGESTED = 0.3
+ ! number of points per minimum wavelength for minimum period estimate
+ real(kind=CUSTOM_REAL),parameter :: NPTS_PER_WAVELENGTH = 5
+
+ !********************************************************************************
+ !real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
+
+
logical :: has_vs_zero
-! initializations
+ ! initializations
if( DT <= 0.0d0) then
DT_PRESENT = .false.
else
@@ -80,6 +88,9 @@
distance_min_glob = HUGEVAL
distance_max_glob = -HUGEVAL
+ elemsize_min_glob = HUGEVAL
+ elemsize_max_glob = -HUGEVAL
+
cmax_glob = -HUGEVAL
pmax_glob = -HUGEVAL
@@ -87,44 +98,13 @@
has_vs_zero = .false.
-! checks courant number & minimum resolved period for each grid cell
+ ! checks courant number & minimum resolved period for each grid cell
do ispec=1,NSPEC_AB
-! determines minimum/maximum velocities within this element
- vpmin = HUGEVAL
- vpmax = -HUGEVAL
- vsmin = HUGEVAL
- vsmax = -HUGEVAL
- ! vp
- where( rho_vp(:,:,:,ispec) > TINYVAL )
- vp_elem(:,:,:) = (FOUR_THIRDS * mustore(:,:,:,ispec) + kappastore(:,:,:,ispec)) / rho_vp(:,:,:,ispec)
- elsewhere
- vp_elem(:,:,:) = 0.0
- endwhere
- ! vs
- where( rho_vs(:,:,:,ispec) > TINYVAL )
- vs_elem(:,:,:) = mustore(:,:,:,ispec) / rho_vs(:,:,:,ispec)
- elsewhere
- vs_elem(:,:,:) = 0.0
- endwhere
+ ! determines minimum/maximum velocities within this element
+ call get_vpvs_minmax(vpmin,vpmax,vsmin,vsmax,ispec,has_vs_zero, &
+ NSPEC_AB,kappastore,mustore,rho_vp,rho_vs)
- val_min = minval(vp_elem(:,:,:))
- val_max = maxval(vp_elem(:,:,:))
-
- vpmin = min(vpmin,val_min(1))
- vpmax = max(vpmax,val_max(1))
-
- val_min = minval(vs_elem(:,:,:))
- val_max = maxval(vs_elem(:,:,:))
-
- ! ignore fluid regions with Vs = 0
- if( val_min(1) > 0.0001 ) then
- vsmin = min(vsmin,val_min(1))
- else
- has_vs_zero = .true.
- endif
- vsmax = max(vsmax,val_max(1))
-
! min/max for whole cpu partition
vpmin_glob = min ( vpmin_glob, vpmin)
vpmax_glob = max ( vpmax_glob, vpmax)
@@ -132,43 +112,23 @@
vsmin_glob = min ( vsmin_glob, vsmin)
vsmax_glob = max ( vsmax_glob, vsmax)
-! compute minimum and maximum distance of GLL points in this grid cell
- distance_min = HUGEVAL
- distance_max = -HUGEVAL
+ ! computes minimum and maximum distance of neighbor GLL points in this grid cell
+ call get_GLL_minmaxdistance(distance_min,distance_max,ispec, &
+ NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
+
+ distance_min_glob = min( distance_min_glob, distance_min)
+ distance_max_glob = max( distance_max_glob, distance_max)
- ! loops over all GLL points
- do k=1,NGLLZ-1
- do j=1,NGLLY-1
- do i=1,NGLLX-1
- iglob_a = ibool(i,j,k,ispec)
+ ! computes minimum and maximum size of this grid cell
+ call get_elem_minmaxsize(elemsize_min,elemsize_max,ispec, &
+ NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
- ! loops over nearest neighbor points
- ! maybe a faster method could be found...
- do kk=k-1,k+1
- do jj=j-1,j+1
- do ii=i-1,i+1
- if( ii < 1 .or. jj < 1 .or. kk < 1 ) cycle
- ! distances between points
- iglob_b = ibool(ii,jj,kk,ispec)
- if( iglob_a /= iglob_b) then
- dx = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
- + ( ystore(iglob_a) - ystore(iglob_b) )**2 &
- + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
- if( dx < distance_min) distance_min = dx
- if( dx > distance_max) distance_max = dx
- endif
- enddo
- enddo
- enddo
+ elemsize_min_glob = min( elemsize_min_glob, elemsize_min)
+ elemsize_max_glob = max( elemsize_max_glob, elemsize_max)
- enddo
- enddo
- enddo
-
- distance_min_glob = min( distance_min_glob, distance_min)
- distance_max_glob = max( distance_max_glob, distance_max)
-
! courant number
+ ! based on minimum GLL point distance and maximum velocity
+ ! i.e. on the maximum ratio of ( velocity / gridsize )
if( DT_PRESENT ) then
cmax = max( vpmax,vsmax ) * DT / distance_min
cmax_glob = max(cmax_glob,cmax)
@@ -179,9 +139,30 @@
dt_suggested_glob = min( dt_suggested_glob, dt_suggested)
! estimation of minimum period resolved
- pmax = distance_max / min( vpmin,vsmin ) * NELEM_PER_WAVELENGTH
+ ! based on average GLL distance within element and minimum velocity
+ !
+ ! rule of thumb (Komatitsch et al. 2005):
+ ! "average number of points per minimum wavelength in an element should be around 5."
+
+ ! average distance between GLL points within this element
+ avg_distance = elemsize_max / NGLLX ! since NGLLX = NGLLY = NGLLZ
+
+ ! biggest possible minimum period such that number of points per minimum wavelength
+ ! npts = ( min(vpmin,vsmin) * pmax ) / avg_distance is about ~ NPTS_PER_WAVELENGTH
+ !
+ ! note: obviously, this estimation depends on the choice of points per wavelength
+ ! which is empirical at the moment.
+ ! also, keep in mind that the minimum period is just an estimation and
+ ! there is no such sharp cut-off period for valid synthetics.
+ ! seismograms become just more and more inaccurate for periods shorter than this estimate.
+ pmax = avg_distance / min( vpmin,vsmin ) * NPTS_PER_WAVELENGTH
pmax_glob = max(pmax_glob,pmax)
+
+ ! old: based on GLL distance, i.e. on maximum ratio ( gridspacing / velocity )
+ !pmax = distance_max / min( vpmin,vsmin ) * NELEM_PER_WAVELENGTH
+ !pmax_glob = max(pmax_glob,pmax)
+
enddo
! determines global min/max values from all cpu partitions
@@ -213,34 +194,46 @@
call min_all_cr(vsmin,vsmin_glob)
call max_all_cr(vsmax,vsmax_glob)
+ ! checks velocities
+ if( vpmin_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vp minimum velocity")
+ endif
+ if( vpmax_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vp maximum velocity")
+ endif
+ if( vsmin_glob < 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vs minimum velocity")
+ endif
+ if( vsmax_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vs maximum velocity")
+ endif
+
! GLL point distance
distance_min = distance_min_glob
distance_max = distance_max_glob
call min_all_cr(distance_min,distance_min_glob)
call max_all_cr(distance_max,distance_max_glob)
+ ! 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)
-! checks mesh
+ ! checks mesh
if( distance_min_glob <= 0.0_CUSTOM_REAL ) then
call exit_mpi(myrank,"error: GLL points minimum distance")
endif
if( distance_max_glob >= HUGEVAL ) then
call exit_mpi(myrank,"error: GLL points maximum distance")
endif
- if( vpmin_glob <= 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vp minimum velocity")
+ if( elemsize_min_glob <= 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: element minimum size")
endif
- if( vpmax_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vp maximum velocity")
- endif
- if( vsmin_glob < 0.0_CUSTOM_REAL ) then
- call exit_mpi(myrank,"error: vs minimum velocity")
- endif
- if( vsmax_glob >= HUGEVAL ) then
- call exit_mpi(myrank,"error: vs maximum velocity")
- endif
+ if( elemsize_max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: element maximum size")
+ endif
-
!! DK DK May 2009: added this to print the minimum and maximum number of elements
!! DK DK May 2009: and points in the CUBIT + SCOTCH mesh
call min_all_i(NSPEC_AB,NSPEC_AB_global_min)
@@ -281,6 +274,9 @@
write(IMAIN,*) '*** Max GLL point distance = ',distance_max_glob
write(IMAIN,*) '*** Min GLL point distance = ',distance_min_glob
write(IMAIN,*) '*** Max/min ratio = ',distance_max_glob/distance_min_glob
+ write(IMAIN,*) '*** Max element size = ',elemsize_max_glob
+ write(IMAIN,*) '*** Min element size = ',elemsize_min_glob
+ write(IMAIN,*) '*** Max/min ratio = ',elemsize_max_glob/elemsize_min_glob
write(IMAIN,*)
write(IMAIN,*) '*** Minimum period resolved = ',pmax_glob
write(IMAIN,*) '*** Maximum suggested time step = ',dt_suggested_glob
@@ -307,5 +303,208 @@
call bcast_all_cr(min_resolved_period,1)
- end subroutine
+ end subroutine check_mesh_resolution
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine get_vpvs_minmax(vpmin,vpmax,vsmin,vsmax,ispec,has_vs_zero, &
+ NSPEC_AB,kappastore,mustore,rho_vp,rho_vs)
+
+! calculates the min/max size of the specified element (ispec)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) :: vpmin,vpmax,vsmin,vsmax
+
+ integer :: ispec
+ logical :: has_vs_zero
+
+ integer :: NSPEC_AB
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ kappastore,mustore,rho_vp,rho_vs
+
+ ! local parameters
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ):: vp_elem,vs_elem
+ real(kind=CUSTOM_REAL), dimension(1) :: val_min,val_max
+
+ ! initializes
+ vpmin = HUGEVAL
+ vpmax = -HUGEVAL
+ vsmin = HUGEVAL
+ vsmax = -HUGEVAL
+
+ ! vp
+ where( rho_vp(:,:,:,ispec) > TINYVAL )
+ vp_elem(:,:,:) = (FOUR_THIRDS * mustore(:,:,:,ispec) &
+ + kappastore(:,:,:,ispec)) / rho_vp(:,:,:,ispec)
+ elsewhere
+ vp_elem(:,:,:) = 0.0
+ endwhere
+
+ val_min = minval(vp_elem(:,:,:))
+ val_max = maxval(vp_elem(:,:,:))
+
+ vpmin = min(vpmin,val_min(1))
+ vpmax = max(vpmax,val_max(1))
+
+ ! vs
+ where( rho_vs(:,:,:,ispec) > TINYVAL )
+ vs_elem(:,:,:) = mustore(:,:,:,ispec) / rho_vs(:,:,:,ispec)
+ elsewhere
+ vs_elem(:,:,:) = 0.0
+ endwhere
+
+ val_min = minval(vs_elem(:,:,:))
+ val_max = maxval(vs_elem(:,:,:))
+
+ ! ignore fluid regions with Vs = 0
+ if( val_min(1) > 0.0001 ) then
+ vsmin = min(vsmin,val_min(1))
+ else
+ has_vs_zero = .true.
+ endif
+ vsmax = max(vsmax,val_max(1))
+
+ end subroutine get_vpvs_minmax
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine get_GLL_minmaxdistance(distance_min,distance_max,ispec, &
+ NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
+
+! calculates the min/max distances between neighboring GLL points within
+! the specified element (ispec)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) :: distance_min,distance_max
+
+ integer :: ispec
+
+ integer :: NSPEC_AB,NGLOB_AB
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: dx,x0,y0,z0
+ integer :: i,j,k,ii,jj,kk,iglob_a,iglob_b
+
+ ! initializes
+ distance_min = HUGEVAL
+ distance_max = -HUGEVAL
+
+ ! loops over all GLL points
+ do k=1,NGLLZ-1
+ do j=1,NGLLY-1
+ do i=1,NGLLX-1
+ iglob_a = ibool(i,j,k,ispec)
+ x0 = xstore(iglob_a)
+ y0 = ystore(iglob_a)
+ z0 = zstore(iglob_a)
+
+ ! loops over nearest neighbor points
+ ! maybe a faster method could be found...
+ do kk=k-1,k+1
+ do jj=j-1,j+1
+ do ii=i-1,i+1
+ if( ii < 1 .or. jj < 1 .or. kk < 1 ) cycle
+ ! distances between points
+ iglob_b = ibool(ii,jj,kk,ispec)
+ if( iglob_a /= iglob_b) then
+ dx = sqrt( ( x0 - xstore(iglob_b) )**2 &
+ + ( y0 - ystore(iglob_b) )**2 &
+ + ( z0 - zstore(iglob_b) )**2 )
+ if( dx < distance_min) distance_min = dx
+ if( dx > distance_max) distance_max = dx
+ endif
+ enddo
+ enddo
+ enddo
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine get_GLL_minmaxdistance
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine get_elem_minmaxsize(elemsize_min,elemsize_max,ispec, &
+ NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
+
+! calculates the min/max size of the specified element (ispec)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max
+
+ integer :: ispec
+
+ integer :: NSPEC_AB,NGLOB_AB
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: dx,x0,y0,z0
+ integer :: i,j,k,icorner,jcorner,iglob_a,iglob_b
+
+ ! corners indices of reference cube faces
+ ! shapes of arrays below
+ integer,dimension(2),parameter :: corner_shape = (/3,NGNOD/)
+ integer,dimension(3,NGNOD),parameter :: corner_ijk = &
+ reshape((/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ, &
+ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),corner_shape)
+
+ ! initializes
+ elemsize_min = HUGEVAL
+ elemsize_max = -HUGEVAL
+
+ ! loops over corners
+ do icorner=1,NGNOD
+ i = corner_ijk(1,icorner)
+ j = corner_ijk(2,icorner)
+ k = corner_ijk(3,icorner)
+
+ iglob_a = ibool(i,j,k,ispec)
+ x0 = xstore(iglob_a)
+ y0 = ystore(iglob_a)
+ z0 = zstore(iglob_a)
+
+ ! loops over all other corners
+ do jcorner = icorner+1,NGNOD
+ i = corner_ijk(1,jcorner)
+ j = corner_ijk(2,jcorner)
+ k = corner_ijk(3,jcorner)
+
+ ! coordinates
+ iglob_b = ibool(i,j,k,ispec)
+
+ ! distances between points
+ if( iglob_a /= iglob_b) then
+ dx = sqrt( ( x0 - xstore(iglob_b) )**2 &
+ + ( y0 - ystore(iglob_b) )**2 &
+ + ( z0 - zstore(iglob_b) )**2 )
+ if( dx < elemsize_min) elemsize_min = dx
+ if( dx > elemsize_max) elemsize_max = dx
+ endif
+
+ enddo
+ enddo
+
+ end subroutine get_elem_minmaxsize
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -69,19 +69,12 @@
logical, dimension(:),allocatable :: mask_ibool_asteroid
- integer :: ixmin, ixmax
- integer :: iymin, iymax
- integer :: izmin, izmax
+ integer :: ixmin, ixmax, iymin, iymax, izmin, izmax
integer, dimension(ngnode) :: n
integer :: e1, e2, e3, e4
- integer :: type
- integer :: ispec
-
- integer :: k
+ integer :: ispec,k,ix,iy,iz,ier,itype,iglob
integer :: npoin_interface_asteroid
- integer :: ix,iy,iz,ier
-
! initializes
allocate( mask_ibool_asteroid(npoin), stat=ier); if( ier /= 0) stop 'error allocating array'
@@ -98,7 +91,7 @@
! spectral element on interface
ispec = my_interfaces(1,ispec_interface,num_interface)
! type of interface: (1) corner point, (2) edge, (4) face
- type = my_interfaces(2,ispec_interface,num_interface)
+ itype = my_interfaces(2,ispec_interface,num_interface)
! gets spectral element corner indices (defines all nodes of face/edge)
do k = 1, ngnode
n(k) = knods(k,ispec)
@@ -111,20 +104,23 @@
e4 = my_interfaces(6,ispec_interface,num_interface)
! gets i,j,k ranges for interface type
- call get_edge(ngnode, n, type, e1, e2, e3, e4, ixmin, ixmax, iymin, iymax, izmin, izmax)
+ call get_edge(ngnode, n, itype, e1, e2, e3, e4, &
+ ixmin, ixmax, iymin, iymax, izmin, izmax)
! counts number and stores indices of (global) points on MPI interface
do iz = min(izmin,izmax), max(izmin,izmax)
do iy = min(iymin,iymax), max(iymin,iymax)
do ix = min(ixmin,ixmax), max(ixmin,ixmax)
+ ! global index
+ iglob = ibool(ix,iy,iz,ispec)
+
! stores global index of point on interface
- if(.not. mask_ibool_asteroid(ibool(ix,iy,iz,ispec))) then
+ if(.not. mask_ibool_asteroid(iglob)) then
! masks point as being accounted for
- mask_ibool_asteroid(ibool(ix,iy,iz,ispec)) = .true.
+ mask_ibool_asteroid(iglob) = .true.
! adds point to interface
npoin_interface_asteroid = npoin_interface_asteroid + 1
- ibool_interfaces_asteroid(npoin_interface_asteroid,num_interface) = &
- ibool(ix,iy,iz,ispec)
+ ibool_interfaces_asteroid(npoin_interface_asteroid,num_interface) = iglob
end if
end do
end do
@@ -145,9 +141,11 @@
!----
!
-subroutine get_edge ( ngnode, n, type, e1, e2, e3, e4, ixmin, ixmax, iymin, iymax, izmin, izmax )
+subroutine get_edge ( ngnode, n, itype, e1, e2, e3, e4, &
+ ixmin, ixmax, iymin, iymax, izmin, izmax )
-! returns range of local (GLL) point indices i,j,k depending on given type for corner point (1), edge (2) or face (4)
+! returns range of local (GLL) point indices i,j,k depending on given type
+! for corner point (1), edge (2) or face (4)
implicit none
@@ -158,7 +156,7 @@
integer, dimension(ngnode), intent(in) :: n
! interface type & nodes
- integer, intent(in) :: type, e1, e2, e3, e4
+ integer, intent(in) :: itype, e1, e2, e3, e4
! local (GLL) i,j,k index ranges
integer, intent(out) :: ixmin, ixmax, iymin, iymax, izmin, izmax
@@ -168,7 +166,7 @@
integer :: valence, i
! determines local indexes for corners/edges/faces
- if ( type == 1 ) then
+ if ( itype == 1 ) then
! corner point
@@ -237,7 +235,7 @@
izmax = NGLLZ
end if
- else if ( type == 2 ) then
+ else if ( itype == 2 ) then
! edges
@@ -402,7 +400,7 @@
end if
end if
- else if (type == 4) then
+ else if (itype == 4) then
! face corners
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/sort_array_coordinates.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/sort_array_coordinates.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/sort_array_coordinates.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -27,7 +27,8 @@
! subroutines to sort MPI buffers to assemble between chunks
- subroutine sort_array_coordinates(npointot,x,y,z,ibool,iglob,loc,ifseg,nglob,ind,ninseg,iwork,work)
+ subroutine sort_array_coordinates(npointot,x,y,z,ibool,iglob,loc,ifseg, &
+ nglob,ind,ninseg,iwork,work)
! this routine MUST be in double precision to avoid sensitivity
! to roundoff errors in the coordinates of the points
@@ -46,7 +47,7 @@
double precision x(npointot),y(npointot),z(npointot)
integer iwork(npointot)
double precision work(npointot)
-
+ ! local parameters
integer ipoin,i,j
integer nseg,ioff,iseg,ig
double precision xtol
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -76,7 +76,8 @@
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
- buffer_send_vector_ext_mesh(:,ipoin,iinterface) = array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
+ buffer_send_vector_ext_mesh(:,ipoin,iinterface) = &
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
enddo
enddo
@@ -105,7 +106,8 @@
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
- array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) &
+ + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
enddo
enddo
@@ -165,7 +167,8 @@
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
- buffer_send_vector_ext_mesh(:,ipoin,iinterface) = array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
+ buffer_send_vector_ext_mesh(:,ipoin,iinterface) = &
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
enddo
enddo
@@ -235,7 +238,8 @@
do iinterface = 1, num_interfaces_ext_mesh
do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
- array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) &
+ + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2011-02-23 21:06:29 UTC (rev 17958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2011-02-23 21:22:42 UTC (rev 17959)
@@ -951,8 +951,16 @@
if( stutm_y >= LATITUDE_MIN .and. stutm_y <= LATITUDE_MAX .and. &
stutm_x >= LONGITUDE_MIN .and. stutm_x <= LONGITUDE_MAX) then
- write(IOUT,*) trim(station_name),' ',trim(network_name),' ',sngl(stlat), &
- ' ',sngl(stlon), ' ',sngl(stele), ' ',sngl(stbur)
+
+ ! w/out formating
+ ! write(IOUT,*) trim(station_name),' ',trim(network_name),' ',sngl(stlat), &
+ ! ' ',sngl(stlon), ' ',sngl(stele), ' ',sngl(stbur)
+
+ ! w/ specific format
+ write(IOUT,'(a10,1x,a10,4e18.6)') &
+ trim(station_name),trim(network_name), &
+ sngl(stlat),sngl(stlon),sngl(stele),sngl(stbur)
+
endif
end if
enddo
More information about the CIG-COMMITS
mailing list