[cig-commits] r8524 - in seismo/2D/SPECFEM2D/trunk: . DATA
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:53:27 PST 2007
Author: walter
Date: 2007-12-07 15:53:26 -0800 (Fri, 07 Dec 2007)
New Revision: 8524
Added:
seismo/2D/SPECFEM2D/trunk/TODO_list
seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/constants_unstruct.h
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/scotchf.h
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/checkgrid.F90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/create_color_image.f90
seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
seismo/2D/SPECFEM2D/trunk/locate_receivers.F90
seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
parallel and unstructured mesher and solver. A TODO list has been included, the readme has been updated (more to come).
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-07 23:53:26 UTC (rev 8524)
@@ -3,6 +3,19 @@
title = Test for M2 UPPA
interfacesfile = interfaces_M2_UPPA_curved.dat
+# data concerning mesh, when generated using third-party app (more info in README)
+read_external_mesh = .false.
+mesh_file = ./DATA/celine_coupe_eros/mesh # file containing the mesh
+nodes_coords_file = ./DATA/celine_coupe_eros/nodes_coords # file containing the nodes coordinates
+materials_file = ./DATA/celine_coupe_eros/mat # file containing the material number for each element
+free_surface_file = ./DATA/celine_coupe_eros/surface_free # file containing the free surface
+absorbing_surface_file = ./DATA/celine_coupe_eros/surface_abs # file containing the absorbing surface
+
+# parameters concerning partitionning
+nproc = 8 # number of processes
+partionning_method = 1 # ascending order = 1, Metis = 2, Scotch = 3
+partitionning_strategy = 01110 #b{strat=m{asc=g,low=g,rat=1.0,type=h,vert=4}} #01110 # options concerning partitionning strategy.
+
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 4000.d0 # abscissa of right side of the model
@@ -14,7 +27,8 @@
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
# absorbing boundaries parameters
-absorbbottom = .true. # absorbing boundary active or not
+absorbing_conditions = .true. # absorbing boundary active or not
+absorbbottom = .true.
absorbright = .true.
absorbtop = .false.
absorbleft = .true.
@@ -71,7 +85,7 @@
# set vs to zero to make a given model acoustic
# the mesh can contain both acoustic and elastic models simultaneously
1 1 2700.d0 3000.d0 1732.051d0 0 0
-2 1 2500.d0 2700.d0 0 0 0 #1558.89d0 0 0
+2 1 2500.d0 2700.d0 0.d0 0 0 #1558.89d0 0 0
3 1 2200.d0 2500.d0 1443.375d0 0 0
4 1 2200.d0 2200.d0 1343.375d0 0 0
# define the different regions of the model in the (nx,nz) spectral element mesh
@@ -79,4 +93,4 @@
1 80 1 20 1
1 80 21 40 2
1 80 41 60 3
-60 70 21 40 4
+60 70 21 40 4
\ No newline at end of file
Modified: seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2007-12-07 23:53:26 UTC (rev 8524)
@@ -1,5 +1,12 @@
- 4
-S001 AA 5000.0000000 15000.0000000 0.0 0.0
-S002 AA 10000.0000000 15000.0000000 0.0 0.0
-S003 AA 15000.0000000 15000.0000000 0.0 0.0
-S004 AA 20000.0000000 15000.0000000 0.0 0.0
+ 11
+S0001 AA 300.0000000 2997.7298909 0.0 0.0
+S0002 AA 640.0000000 3008.0430011 0.0 0.0
+S0003 AA 980.0000000 3090.8224062 0.0 0.0
+S0004 AA 1320.0000000 3283.0303923 0.0 0.0
+S0005 AA 1660.0000000 3347.8768862 0.0 0.0
+S0006 AA 2000.0000000 3250.0000000 0.0 0.0
+S0007 AA 2340.0000000 3197.3138031 0.0 0.0
+S0008 AA 2680.0000000 3150.9619873 0.0 0.0
+S0009 AA 3020.0000000 3086.5939051 0.0 0.0
+S0010 AA 3360.0000000 3042.8523748 0.0 0.0
+S0011 AA 3700.0000000 3020.6886768 0.0 0.0
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:53:26 UTC (rev 8524)
@@ -18,22 +18,26 @@
#FLAGS_CHECK = $(FLAGS_NOCHECK) -check bounds
# GNU gfortran
-F90 = gfortran
+F90 = /opt/openmpi-1.2.1/gfortran64/bin/mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
+#F90 = gfortran
#FLAGS_NOCHECK = -O3 -march=opteron -m64 -mfpmath=sse,387
FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O2 -Wunused-labels -Waliasing -Wampersand -Wsurprising -Wline-truncation -Wunderflow
FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
LINK = $(F90)
-OBJS_MESHFEM2D = $O/meshfem2D.o $O/read_value_parameters.o
+LIB = /opt/metis-4.0.1/gcc64/lib/libmetis.a /opt/scotch-4.0/gcc64/lib/libscotch.a /opt/scotch-4.0/gcc64/lib/libscotcherr.a
+OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o
+
OBJS_SPECFEM2D = $O/checkgrid.o $O/datim.o $O/enforce_acoustic_free_surface.o\
$O/compute_forces_acoustic.o $O/compute_forces_elastic.o\
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivation_matrices.o\
$O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o\
$O/specfem2D.o $O/write_seismograms.o $O/define_external_model.o $O/createnum_fast.o $O/createnum_slow.o\
$O/define_shape_functions.o $O/create_color_image.o $O/compute_vector_field.o $O/compute_pressure.o\
- $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/numerical_recipes.o
+ $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/numerical_recipes.o\
+ $O/construct_acoustic_surface.o $O/assemble_MPI.o
default: clean meshfem2D specfem2D convolve_source_timefunction
@@ -43,7 +47,7 @@
/bin/rm -r -f xmeshfem2D xmeshfem2D.trace xspecfem2D xspecfem2D.trace $O/*.o *.o $O/*.il *.mod core xconvolve_source_timefunction
meshfem2D: $(OBJS_MESHFEM2D)
- $(LINK) $(FLAGS_CHECK) -o xmeshfem2D $(OBJS_MESHFEM2D)
+ $(LINK) $(FLAGS_CHECK) -o xmeshfem2D $(OBJS_MESHFEM2D) $(LIB)
### use optimized compilation option for solver only
specfem2D: $(OBJS_SPECFEM2D)
@@ -147,3 +151,11 @@
$O/write_seismograms.o: write_seismograms.F90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/write_seismograms.o write_seismograms.F90
+$O/part_unstruct.o: part_unstruct.F90 constants_unstruct.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/part_unstruct.o part_unstruct.F90
+
+$O/construct_acoustic_surface.o: construct_acoustic_surface.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/construct_acoustic_surface.o construct_acoustic_surface.f90
+
+$O/assemble_MPI.o: assemble_MPI.F90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/assemble_MPI.o assemble_MPI.F90
Added: seismo/2D/SPECFEM2D/trunk/TODO_list
===================================================================
--- seismo/2D/SPECFEM2D/trunk/TODO_list 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/TODO_list 2007-12-07 23:53:26 UTC (rev 8524)
@@ -0,0 +1,16 @@
+- enabling output of vectxxx.ps files in parallel.
+- enabling output of source.txt file in parallel.
+- looking for a better way to show mesh partitioning.
+- splitting file part_unstruct.F90 in several files for clarity purpose.
+- improving compiling with SCOTCH (issue with header file scotchf.h which is Fortran77 legal). Having our own scotchf.h file (without the comments) is ugly.
+- comparing the different partitioning methods for METIS and SCOTCH, and finding a good default for SCOTCH.
+- partitioning using weights for load-balancing.
+- choosing a way to use assign_external_model?
+- adding comments.
+- adding option to create or not STATIONS file (useful for complex geometries like asteroids).
+- checking for points with different normals for absorbing conditions, when the absorbing edges are not in the same elements (similar to what is done for the corners).
+- scripts for translating GID/CUBIT meshes into files for xmeshfem2D.
+- modifying scripts for UPPA cluster (when FS sync issues are solved and remote commands are available).
+- manual for unstructured meshes.
+- getting rid of constants_unstruct.h.
+- enabling use of real or double precision.
\ No newline at end of file
Added: seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -0,0 +1,760 @@
+subroutine prepare_assemble_MPI (myrank,nspec,ibool, &
+ knods, ngnod, &
+ npoin, elastic, &
+ ninterface, max_interface_size, &
+ my_neighbours, my_nelmnts_neighbours, my_interfaces, &
+ ibool_interfaces_acoustic, ibool_interfaces_elastic, &
+ nibool_interfaces_acoustic, nibool_interfaces_elastic, &
+ inum_interfaces_acoustic, inum_interfaces_elastic, &
+ ninterface_acoustic, ninterface_elastic &
+ )
+
+
+ implicit none
+ include 'constants.h'
+
+
+ integer, intent(in) :: nspec, myrank, npoin, ngnod
+ logical, dimension(nspec), intent(in) :: elastic
+ integer, dimension(ngnod,nspec), intent(in) :: knods
+ integer, dimension(NGLLX,NGLLZ,nspec), intent(in) :: ibool
+
+ !integer, dimension(nspec) :: inner_to_glob_ispec
+ !integer, dimension(nspec) :: interface_to_glob_ispec
+
+ !integer, intent(inout) :: nspec_inner_known
+ !integer, intent(inout) :: nspec_interface_known
+
+
+ integer :: ninterface
+ integer :: max_interface_size
+ integer, dimension(ninterface) :: my_neighbours
+ integer, dimension(ninterface) :: my_nelmnts_neighbours
+ integer, dimension(4,max_interface_size,ninterface) :: my_interfaces
+ integer, dimension(NGLLX*max_interface_size,ninterface) :: &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic
+ integer, dimension(ninterface) :: &
+ nibool_interfaces_acoustic,nibool_interfaces_elastic
+
+
+ integer, dimension(ninterface), intent(out) :: &
+ inum_interfaces_acoustic, inum_interfaces_elastic
+ integer, intent(out) :: ninterface_acoustic, ninterface_elastic
+
+ integer :: num_interface
+ integer :: ispec_interface
+
+ logical, dimension(nspec) :: ispec_is_inner
+ logical, dimension(npoin) :: mask_ibool_acoustic
+ logical, dimension(npoin) :: mask_ibool_elastic
+
+ integer :: ixmin, ixmax
+ integer :: izmin, izmax
+ integer, dimension(ngnod) :: n
+ integer :: e1, e2
+ integer :: type
+ integer :: ispec
+
+ integer :: k
+ integer :: npoin_interface_acoustic
+ integer :: npoin_interface_elastic
+
+ integer :: ix,iz
+ integer :: ier
+
+ integer :: sens
+
+
+ ibool_interfaces_acoustic(:,:) = 0
+ nibool_interfaces_acoustic(:) = 0
+ ibool_interfaces_elastic(:,:) = 0
+ nibool_interfaces_elastic(:) = 0
+
+ do num_interface = 1, ninterface
+ npoin_interface_acoustic = 0
+ npoin_interface_elastic = 0
+ mask_ibool_acoustic(:) = .false.
+ mask_ibool_elastic(:) = .false.
+
+ do ispec_interface = 1, my_nelmnts_neighbours(num_interface)
+ ispec = my_interfaces(1,ispec_interface,num_interface)
+ type = my_interfaces(2,ispec_interface,num_interface)
+ do k = 1, ngnod
+ n(k) = knods(k,ispec)
+ end do
+ e1 = my_interfaces(3,ispec_interface,num_interface)
+ e2 = my_interfaces(4,ispec_interface,num_interface)
+
+ call get_edge(ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax, sens)
+
+ do iz = izmin, izmax, sens
+ do ix = ixmin, ixmax, sens
+
+ if ( elastic(ispec) ) then
+
+ if(.not. mask_ibool_elastic(ibool(ix,iz,ispec))) then
+ mask_ibool_elastic(ibool(ix,iz,ispec)) = .true.
+ npoin_interface_elastic = npoin_interface_elastic + 1
+ ibool_interfaces_elastic(npoin_interface_elastic,num_interface)=&
+ ibool(ix,iz,ispec)
+ end if
+ else
+ if(.not. mask_ibool_acoustic(ibool(ix,iz,ispec))) then
+ mask_ibool_acoustic(ibool(ix,iz,ispec)) = .true.
+ npoin_interface_acoustic = npoin_interface_acoustic + 1
+ ibool_interfaces_acoustic(npoin_interface_acoustic,num_interface)=&
+ ibool(ix,iz,ispec)
+ end if
+ end if
+ end do
+ end do
+
+ end do
+ nibool_interfaces_acoustic(num_interface) = npoin_interface_acoustic
+ nibool_interfaces_elastic(num_interface) = npoin_interface_elastic
+
+ end do
+
+
+ ninterface_acoustic = 0
+ ninterface_elastic = 0
+ do num_interface = 1, ninterface
+ if ( nibool_interfaces_acoustic(num_interface) > 0 ) then
+ ninterface_acoustic = ninterface_acoustic + 1
+ inum_interfaces_acoustic(ninterface_acoustic) = num_interface
+ end if
+ if ( nibool_interfaces_elastic(num_interface) > 0 ) then
+ ninterface_elastic = ninterface_elastic + 1
+ inum_interfaces_elastic(ninterface_elastic) = num_interface
+ end if
+ end do
+
+
+!!$ nspec_inner_known = 0
+!!$ do ispec = 1, nspec
+!!$ if ( ispec_is_inner(ispec) ) then
+!!$ nspec_inner_known = nspec_inner_known + 1
+!!$ end if
+!!$ end do
+
+ !allocate(inner_to_glob_ispec(nspec_inner_known))
+ !allocate(interface_to_glob_ispec(nspec-nspec_inner_known))
+
+
+
+!!$ nspec_inner_known = 0
+!!$ nspec_interface_known = 0
+!!$ do ispec = 1, nspec
+!!$ if ( ispec_is_inner(ispec) ) then
+!!$ nspec_inner_known = nspec_inner_known + 1
+!!$ inner_to_glob_ispec(nspec_inner_known) = ispec
+!!$ else
+!!$ nspec_interface_known = nspec_interface_known + 1
+!!$ interface_to_glob_ispec(nspec_interface_known) = ispec
+!!$ end if
+!!$
+!!$ end do
+
+
+end subroutine prepare_assemble_MPI
+
+
+
+subroutine get_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax, sens )
+
+ implicit none
+ include "constants.h"
+
+ integer, intent(in) :: ngnod
+ integer, dimension(ngnod), intent(in) :: n
+ integer, intent(in) :: type, e1, e2
+ integer, intent(out) :: ixmin, ixmax, izmin, izmax
+ integer, intent(out) :: sens
+
+ if ( type == 1 ) then
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ ixmax = 1
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ ixmax = 1
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ sens = 1
+ else
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ izmin = 1
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ izmax = 1
+ sens = 1
+ end if
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ izmax = NGLLZ
+ sens = 1
+ end if
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ izmin = 1
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ izmax = NGLLZ
+ sens = 1
+ end if
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ izmax = 1
+ sens = -1
+ end if
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ izmin = NGLLZ
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ izmax = NGLLZ
+ sens = -1
+ end if
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ izmax = 1
+ sens = -1
+ end if
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ izmin = NGLLZ
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ izmax = 1
+ sens = -1
+ end if
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ izmax = NGLLZ
+ sens = 1
+ end if
+ end if
+ end if
+
+
+end subroutine get_edge
+
+
+
+#ifdef USE_MPI
+
+subroutine create_MPI_requests_SEND_RECV_acoustic(myrank, &
+ ninterface, ninterface_acoustic, &
+ nibool_interfaces_acoustic, &
+ my_neighbours, &
+ max_ibool_interfaces_size_acoustic, &
+ buffer_send_faces_vector_acoustic, &
+ buffer_recv_faces_vector_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ inum_interfaces_acoustic &
+ )
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: ninterface, ninterface_acoustic
+ integer, dimension(ninterface), intent(in) :: inum_interfaces_acoustic
+ integer, intent(in) :: max_ibool_interfaces_size_acoustic
+ double precision, dimension(max_ibool_interfaces_size_acoustic,ninterface_acoustic), intent(in) :: &
+ buffer_send_faces_vector_acoustic
+ double precision, dimension(max_ibool_interfaces_size_acoustic,ninterface_acoustic), intent(in) :: &
+ buffer_recv_faces_vector_acoustic
+ integer, dimension(ninterface_acoustic*2), intent(inout) :: tab_requests_send_recv_acoustic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic
+ integer, dimension(ninterface), intent(in) :: my_neighbours
+
+
+ integer :: inum_interface,num_interface
+ integer :: ier
+
+
+ do inum_interface = 1, ninterface_acoustic
+
+ num_interface = inum_interfaces_acoustic(inum_interface)
+
+ call MPI_Send_init ( buffer_send_faces_vector_acoustic(1,inum_interface), &
+ nibool_interfaces_acoustic(num_interface), MPI_DOUBLE_PRECISION, &
+ my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
+ tab_requests_send_recv_acoustic(inum_interface), ier)
+ call MPI_Recv_init ( buffer_recv_faces_vector_acoustic(1,inum_interface), &
+ nibool_interfaces_acoustic(num_interface), MPI_DOUBLE_PRECISION, &
+ my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
+ tab_requests_send_recv_acoustic(ninterface_acoustic+inum_interface), ier)
+ end do
+
+
+end subroutine create_MPI_requests_SEND_RECV_acoustic
+
+
+
+subroutine create_MPI_requests_SEND_RECV_elastic(myrank, &
+ ninterface, ninterface_elastic, &
+ nibool_interfaces_elastic, &
+ my_neighbours, &
+ max_ibool_interfaces_size_elastic, &
+ buffer_send_faces_vector_elastic, &
+ buffer_recv_faces_vector_elastic, &
+ tab_requests_send_recv_elastic, &
+ inum_interfaces_elastic &
+ )
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: ninterface, ninterface_elastic
+ integer, dimension(ninterface), intent(in) :: inum_interfaces_elastic
+ integer, intent(in) :: max_ibool_interfaces_size_elastic
+ double precision, dimension(max_ibool_interfaces_size_elastic,ninterface_elastic), intent(in) :: &
+ buffer_send_faces_vector_elastic
+ double precision, dimension(max_ibool_interfaces_size_elastic,ninterface_elastic), intent(in) :: &
+ buffer_recv_faces_vector_elastic
+ integer, dimension(ninterface_elastic*2), intent(inout) :: tab_requests_send_recv_elastic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_elastic
+ integer, dimension(ninterface), intent(in) :: my_neighbours
+
+ integer :: inum_interface,num_interface
+ integer :: ier
+
+
+
+ do inum_interface = 1, ninterface_elastic
+
+ num_interface = inum_interfaces_elastic(inum_interface)
+
+ call MPI_Send_init ( buffer_send_faces_vector_elastic(1,inum_interface), &
+ NDIM*nibool_interfaces_elastic(num_interface), MPI_DOUBLE_PRECISION, &
+ my_neighbours(num_interface), 13, MPI_COMM_WORLD, &
+ tab_requests_send_recv_elastic(inum_interface), ier)
+ call MPI_Recv_init ( buffer_recv_faces_vector_elastic(1,inum_interface), &
+ NDIM*nibool_interfaces_elastic(num_interface), MPI_DOUBLE_PRECISION, &
+ my_neighbours(num_interface), 13, MPI_COMM_WORLD, &
+ tab_requests_send_recv_elastic(ninterface_elastic+inum_interface), ier)
+ end do
+
+
+end subroutine create_MPI_requests_SEND_RECV_elastic
+
+
+
+subroutine assemble_MPI_scalar(myrank,array_val1, array_val2,npoin, &
+ ninterface, max_interface_size, max_ibool_interfaces_size_acoustic, max_ibool_interfaces_size_elastic, &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic, nibool_interfaces_acoustic,nibool_interfaces_elastic, my_neighbours)
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ ! array to assemble
+ double precision, dimension(npoin), intent(inout) :: array_val1, array_val2
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: npoin
+ integer, intent(in) :: ninterface
+ integer, intent(in) :: max_interface_size
+ integer, intent(in) :: max_ibool_interfaces_size_acoustic, max_ibool_interfaces_size_elastic
+ integer, dimension(NGLLX*max_interface_size,ninterface), intent(in) :: &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic,nibool_interfaces_elastic
+ integer, dimension(ninterface), intent(in) :: my_neighbours
+
+
+ double precision, dimension(max_ibool_interfaces_size_acoustic+max_ibool_interfaces_size_elastic, ninterface) :: &
+ buffer_send_faces_scalar, &
+ buffer_recv_faces_scalar
+ integer :: msg_status(MPI_STATUS_SIZE)
+ integer, dimension(ninterface) :: msg_requests
+ integer :: ipoin, num_interface
+ integer :: ier
+
+ integer :: i
+
+
+ do num_interface = 1, ninterface
+
+ print *, 'QQQQQ', myrank,num_interface,nibool_interfaces_acoustic(num_interface), nibool_interfaces_elastic(num_interface)
+
+ ipoin = 0
+ do i = 1, nibool_interfaces_acoustic(num_interface)
+ ipoin = ipoin + 1
+ buffer_send_faces_scalar(ipoin,num_interface) = &
+ array_val1(ibool_interfaces_acoustic(i,num_interface))
+ end do
+
+ do i = 1, nibool_interfaces_elastic(num_interface)
+ ipoin = ipoin + 1
+ buffer_send_faces_scalar(ipoin,num_interface) = &
+ array_val2(ibool_interfaces_elastic(i,num_interface))
+ end do
+
+ call MPI_isend ( buffer_send_faces_scalar(1,num_interface), &
+ nibool_interfaces_acoustic(num_interface)+nibool_interfaces_elastic(num_interface), MPI_DOUBLE_PRECISION, &
+ my_neighbours(num_interface), 11, &
+ MPI_COMM_WORLD, msg_requests(num_interface), ier)
+
+ end do
+
+
+ do num_interface = 1, ninterface
+ call MPI_recv ( buffer_recv_faces_scalar(1,num_interface), &
+ nibool_interfaces_acoustic(num_interface)+nibool_interfaces_elastic(num_interface), MPI_DOUBLE_PRECISION, &
+ my_neighbours(num_interface), 11, &
+ MPI_COMM_WORLD, msg_status(1), ier)
+
+ ipoin = 0
+ do i = 1, nibool_interfaces_acoustic(num_interface)
+ ipoin = ipoin + 1
+ array_val1(ibool_interfaces_acoustic(i,num_interface)) = array_val1(ibool_interfaces_acoustic(i,num_interface)) + &
+ buffer_recv_faces_scalar(ipoin,num_interface)
+ end do
+
+ do i = 1, nibool_interfaces_elastic(num_interface)
+ ipoin = ipoin + 1
+ array_val2(ibool_interfaces_elastic(i,num_interface)) = array_val2(ibool_interfaces_elastic(i,num_interface)) + &
+ buffer_recv_faces_scalar(ipoin,num_interface)
+ end do
+
+ end do
+
+ call MPI_BARRIER(mpi_comm_world,ier)
+
+
+end subroutine assemble_MPI_scalar
+
+
+
+subroutine assemble_MPI_vector_acoustic_start(myrank,array_val1,npoin, &
+ ninterface, ninterface_acoustic, &
+ inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_acoustic,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ buffer_send_faces_vector_acoustic, &
+ buffer_recv_faces_vector_acoustic &
+ )
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ ! array to assemble
+ double precision, dimension(npoin), intent(in) :: array_val1
+
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: npoin
+ integer, intent(in) :: ninterface, ninterface_acoustic
+ integer, dimension(ninterface), intent(in) :: inum_interfaces_acoustic
+ integer, intent(in) :: max_interface_size
+ integer, intent(in) :: max_ibool_interfaces_size_acoustic
+ integer, dimension(NGLLX*max_interface_size,ninterface), intent(in) :: ibool_interfaces_acoustic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic
+ integer, dimension(ninterface_acoustic*2), intent(inout) :: tab_requests_send_recv_acoustic
+ double precision, dimension(max_ibool_interfaces_size_acoustic,ninterface_acoustic), intent(inout) :: &
+ buffer_send_faces_vector_acoustic, buffer_recv_faces_vector_acoustic
+
+ integer :: ipoin, num_interface, inum_interface
+ integer :: ier
+ integer, dimension(MPI_STATUS_SIZE,ninterface_acoustic*2) :: tab_statuses_acoustic
+
+ integer :: i
+
+
+ do inum_interface = 1, ninterface_acoustic
+
+ num_interface = inum_interfaces_acoustic(inum_interface)
+
+ ipoin = 0
+ do i = 1, nibool_interfaces_acoustic(num_interface)
+ ipoin = ipoin + 1
+ buffer_send_faces_vector_acoustic(ipoin,inum_interface) = &
+ array_val1(ibool_interfaces_acoustic(i,num_interface))
+ end do
+
+ end do
+
+ do inum_interface = 1, ninterface_acoustic*2
+ call MPI_START(tab_requests_send_recv_acoustic(inum_interface), ier)
+ if ( ier /= MPI_SUCCESS ) then
+ call exit_mpi(myrank,'MPI_start unsuccessful in assemble_MPI_vector_start')
+ end if
+ end do
+
+!call MPI_Startall ( ninterface*2, tab_requests_send_recv(1), ier )
+
+
+end subroutine assemble_MPI_vector_acoustic_start
+
+
+
+subroutine assemble_MPI_vector_elastic_start(myrank,array_val2,npoin, &
+ ninterface, ninterface_elastic, &
+ inum_interfaces_elastic, &
+ max_interface_size, max_ibool_interfaces_size_elastic,&
+ ibool_interfaces_elastic, nibool_interfaces_elastic, &
+ tab_requests_send_recv_elastic, &
+ buffer_send_faces_vector_elastic, &
+ buffer_recv_faces_vector_elastic &
+ )
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ ! array to assemble
+ double precision, dimension(NDIM,npoin), intent(in) :: array_val2
+
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: npoin
+ integer, intent(in) :: ninterface, ninterface_elastic
+ integer, dimension(ninterface), intent(in) :: inum_interfaces_elastic
+ integer, intent(in) :: max_interface_size
+ integer, intent(in) :: max_ibool_interfaces_size_elastic
+ integer, dimension(NGLLX*max_interface_size,ninterface), intent(in) :: ibool_interfaces_elastic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_elastic
+ integer, dimension(ninterface_elastic*2), intent(inout) :: tab_requests_send_recv_elastic
+ double precision, dimension(max_ibool_interfaces_size_elastic,ninterface_elastic), intent(inout) :: &
+ buffer_send_faces_vector_elastic, buffer_recv_faces_vector_elastic
+
+
+
+ integer :: ipoin, num_interface, inum_interface
+ integer :: ier
+ integer, dimension(MPI_STATUS_SIZE,ninterface_elastic) :: tab_statuses_elastic
+
+ integer :: i
+
+
+ do inum_interface = 1, ninterface_elastic
+
+ num_interface = inum_interfaces_elastic(inum_interface)
+
+ ipoin = 0
+ do i = 1, nibool_interfaces_elastic(num_interface)
+ buffer_send_faces_vector_elastic(ipoin+1:ipoin+2,inum_interface) = &
+ array_val2(:,ibool_interfaces_elastic(i,num_interface))
+ ipoin = ipoin + 2
+ end do
+
+
+ end do
+
+ do inum_interface = 1, ninterface_elastic*2
+ call MPI_START(tab_requests_send_recv_elastic(inum_interface), ier)
+ if ( ier /= MPI_SUCCESS ) then
+ call exit_mpi(myrank,'MPI_start unsuccessful in assemble_MPI_vector_start')
+ end if
+ end do
+
+!call MPI_Startall ( ninterface*2, tab_requests_send_recv(1), ier )
+
+
+end subroutine assemble_MPI_vector_elastic_start
+
+
+
+subroutine assemble_MPI_vector_acoustic_wait(myrank,array_val1,npoin, &
+ ninterface, ninterface_acoustic, &
+ inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_acoustic,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ buffer_send_faces_vector_acoustic, &
+ buffer_recv_faces_vector_acoustic &
+ )
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ ! array to assemble
+ double precision, dimension(npoin), intent(inout) :: array_val1
+
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: npoin
+ integer, intent(in) :: ninterface, ninterface_acoustic
+ integer, dimension(ninterface), intent(in) :: inum_interfaces_acoustic
+ integer, intent(in) :: max_interface_size
+ integer, intent(in) :: max_ibool_interfaces_size_acoustic
+ integer, dimension(NGLLX*max_interface_size,ninterface), intent(in) :: ibool_interfaces_acoustic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic
+ integer, dimension(ninterface_acoustic*2), intent(inout) :: tab_requests_send_recv_acoustic
+ double precision, dimension(max_ibool_interfaces_size_acoustic,ninterface_acoustic), intent(inout) :: &
+ buffer_send_faces_vector_acoustic, buffer_recv_faces_vector_acoustic
+
+ integer :: ipoin, num_interface, inum_interface
+ integer :: ier
+ integer, dimension(MPI_STATUS_SIZE,ninterface_acoustic*2) :: tab_statuses_acoustic
+
+ integer :: i
+
+
+ call MPI_Waitall ( ninterface_acoustic*2, tab_requests_send_recv_acoustic(1), tab_statuses_acoustic(1,1), ier )
+ if ( ier /= MPI_SUCCESS ) then
+ call exit_mpi(myrank,'MPI_WAITALL unsuccessful in assemble_MPI_vector_wait')
+ end if
+
+ do inum_interface = 1, ninterface_acoustic
+
+ num_interface = inum_interfaces_acoustic(inum_interface)
+
+ ipoin = 0
+ do i = 1, nibool_interfaces_acoustic(num_interface)
+ ipoin = ipoin + 1
+ array_val1(ibool_interfaces_acoustic(i,num_interface)) = array_val1(ibool_interfaces_acoustic(i,num_interface)) + &
+ buffer_recv_faces_vector_acoustic(ipoin,inum_interface)
+ end do
+
+ end do
+
+
+end subroutine assemble_MPI_vector_acoustic_wait
+
+
+
+subroutine assemble_MPI_vector_elastic_wait(myrank,array_val2,npoin, &
+ ninterface, ninterface_elastic, &
+ inum_interfaces_elastic, &
+ max_interface_size, max_ibool_interfaces_size_elastic,&
+ ibool_interfaces_elastic, nibool_interfaces_elastic, &
+ tab_requests_send_recv_elastic, &
+ buffer_send_faces_vector_elastic, &
+ buffer_recv_faces_vector_elastic &
+ )
+
+ implicit none
+
+ include 'constants.h'
+ include 'mpif.h'
+
+
+ ! array to assemble
+ double precision, dimension(NDIM,npoin), intent(inout) :: array_val2
+
+
+ integer, intent(in) :: myrank
+ integer, intent(in) :: npoin
+ integer, intent(in) :: ninterface, ninterface_elastic
+ integer, dimension(ninterface), intent(in) :: inum_interfaces_elastic
+ integer, intent(in) :: max_interface_size
+ integer, intent(in) :: max_ibool_interfaces_size_elastic
+ integer, dimension(NGLLX*max_interface_size,ninterface), intent(in) :: ibool_interfaces_elastic
+ integer, dimension(ninterface), intent(in) :: nibool_interfaces_elastic
+ integer, dimension(ninterface_elastic*2), intent(inout) :: tab_requests_send_recv_elastic
+ double precision, dimension(max_ibool_interfaces_size_elastic,ninterface_elastic), intent(inout) :: &
+ buffer_send_faces_vector_elastic, buffer_recv_faces_vector_elastic
+
+ integer :: ipoin, num_interface, inum_interface
+ integer :: ier
+ integer, dimension(MPI_STATUS_SIZE,ninterface_elastic*2) :: tab_statuses_elastic
+
+ integer :: i
+
+
+ call MPI_Waitall ( ninterface_elastic*2, tab_requests_send_recv_elastic(1), tab_statuses_elastic(1,1), ier )
+ if ( ier /= MPI_SUCCESS ) then
+ call exit_mpi(myrank,'MPI_WAITALL unsuccessful in assemble_MPI_vector_wait')
+ end if
+
+ do inum_interface = 1, ninterface_elastic
+
+ num_interface = inum_interfaces_elastic(inum_interface)
+
+ ipoin = 0
+ do i = 1, nibool_interfaces_elastic(num_interface)
+ array_val2(:,ibool_interfaces_elastic(i,num_interface)) = array_val2(:,ibool_interfaces_elastic(i,num_interface)) + &
+ buffer_recv_faces_vector_elastic(ipoin+1:ipoin+2,inum_interface)
+ ipoin = ipoin + 2
+ end do
+
+ end do
+
+
+end subroutine assemble_MPI_vector_elastic_wait
+
+
+
+subroutine exit_MPI(myrank,error_msg)
+
+ implicit none
+
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+
+ ! identifier for error message file
+ integer, parameter :: IERROR = 30
+
+ integer myrank
+ character(len=*) error_msg
+
+ integer ier
+ character(len=80) outputname
+
+ ! write error message to screen
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI... proc ',myrank
+
+ ! write error message to file
+ write(outputname,"('error_message',i6.6,'.txt')") myrank
+ open(unit=IERROR,file=outputname,status='unknown')
+ write(IERROR,*) error_msg(1:len(error_msg))
+ write(IERROR,*) 'Error detected, aborting MPI... proc ',myrank
+ close(IERROR)
+
+
+ ! stop all the MPI processes, and exit
+ ! on some machines, MPI_FINALIZE needs to be called before MPI_ABORT
+ call MPI_FINALIZE(ier)
+ call MPI_ABORT(ier)
+ stop 'error, program ended in exit_MPI'
+
+end subroutine exit_MPI
+
+
+#endif
+
+
+
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -13,13 +13,16 @@
subroutine checkgrid(vpext,vsext,rhoext,density,elastcoef,ibool,kmato,coord,npoin,vpmin,vpmax, &
assign_external_model,nspec,numat,deltat,f0,t0,initialfield,time_function_type, &
- coorg,xinterp,zinterp,shapeint,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic)
+ coorg,xinterp,zinterp,shapeint,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic,myrank,nproc)
! check the mesh, stability and number of points per wavelength
implicit none
include "constants.h"
+#ifdef USE_MPI
+ include 'mpif.h'
+#endif
integer i,j,ispec,material,npoin,nspec,numat,time_function_type
@@ -50,6 +53,28 @@
double precision :: xmax,zmax,height,usoffset,sizex,sizez,courant_stability_number
double precision :: x1,z1,x2,z2,ratio_page,xmin,zmin,lambdaS_local,lambdaP_local
+ double precision :: vpmin_glob,vpmax_glob,vsmin_glob,vsmax_glob,densmin_glob,densmax_glob
+ double precision :: distance_min_glob,distance_max_glob
+ double precision :: courant_stability_number_max_glob,lambdaPmin_glob,lambdaPmax_glob,lambdaSmin_glob,lambdaSmax_glob
+ double precision :: xmin_glob, xmax_glob, zmin_glob, zmax_glob
+ logical :: any_elastic_glob
+ double precision, dimension(2,nspec*5) :: coorg_send
+ double precision, dimension(:,:), allocatable :: coorg_recv
+ integer, dimension(nspec) :: RGB_send
+ integer, dimension(:), allocatable :: RGB_recv
+ real, dimension(nspec) :: greyscale_send
+ real, dimension(:), allocatable :: greyscale_recv
+ integer :: nspec_recv
+ integer :: num_ispec
+ integer :: myrank, iproc, nproc
+ integer :: ier
+ integer :: nproc_per_cat, color_cat
+ real :: RGB_step, R_val, G_val, B_val
+
+#ifdef USE_MPI
+ integer, dimension(MPI_STATUS_SIZE) :: request_mpi_status
+#endif
+
integer knods(ngnod,nspec)
double precision xinterp(pointsdisp,pointsdisp),zinterp(pointsdisp,pointsdisp)
@@ -179,6 +204,39 @@
enddo
+ any_elastic_glob = any_elastic
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (vpmin, vpmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vpmax, vpmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vsmin, vsmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vsmax, vsmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (densmin, densmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (densmax, densmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (distance_min, distance_min_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (distance_max, distance_max_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (courant_stability_number_max, courant_stability_number_max_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaPmin, lambdaPmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaPmax, lambdaPmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaSmin, lambdaSmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaSmax, lambdaSmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (any_elastic, any_elastic_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
+ vpmin = vpmin_glob
+ vpmax = vpmax_glob
+ vsmin = vsmin_glob
+ vsmax = vsmax_glob
+ densmin = densmin_glob
+ densmax = densmax_glob
+ distance_min = distance_min_glob
+ distance_max = distance_max_glob
+ courant_stability_number_max = courant_stability_number_max_glob
+ lambdaPmin = lambdaPmin_glob
+ lambdaPmax = lambdaPmax_glob
+ lambdaSmin = lambdaSmin_glob
+ lambdaSmax = lambdaSmax_glob
+
+#endif
+
write(IOUT,*)
write(IOUT,*) '********'
write(IOUT,*) 'Model: P velocity min,max = ',vpmin,vpmax
@@ -245,12 +303,25 @@
xmax = maxval(coord(1,:))
zmax = maxval(coord(2,:))
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (xmin, xmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (xmax, xmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (zmin, zmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (zmax, zmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ xmin = xmin_glob
+ xmax = xmax_glob
+ zmin = zmin_glob
+ zmax = zmax_glob
+
+#endif
+
! ratio of physical page size/size of the domain meshed
ratio_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin)) / 100.d0
print *
print *,'Creating PostScript file with stability condition'
+ if ( myrank == 0 ) then
!
!---- open PostScript file
!
@@ -354,10 +425,16 @@
write(24,*) '%'
write(24,*) '0 setgray'
+ num_ispec = 0
+ end if
+
do ispec = 1, nspec
+
+ if ( myrank == 0 ) then
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',ispec
+ end if
- write(24,*) '% elem ',ispec
-
do i=1,pointsdisp
do j=1,pointsdisp
xinterp(i,j) = 0.d0
@@ -376,8 +453,13 @@
z1 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x1 = x1 * centim
z1 = z1 * centim
- write(24,*) 'mark'
- write(24,681) x1,z1
+ if ( myrank == 0 ) then
+ write(24,*) 'mark'
+ write(24,681) x1,z1
+ else
+ coorg_send(1,(ispec-1)*5+1) = x1
+ coorg_send(2,(ispec-1)*5+1) = z1
+ end if
! draw straight lines if elements have 4 nodes
@@ -386,7 +468,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+2) = x2
+ coorg_send(2,(ispec-1)*5+2) = z2
+ end if
ir=pointsdisp
is=pointsdisp
@@ -394,7 +481,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+3) = x2
+ coorg_send(2,(ispec-1)*5+3) = z2
+ end if
is=pointsdisp
ir=1
@@ -402,7 +494,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+4) = x2
+ coorg_send(2,(ispec-1)*5+4) = z2
+ end if
ir=1
is=2
@@ -410,10 +507,14 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ write(24,*) 'CO'
+ else
+ coorg_send(1,(ispec-1)*5+5) = x2
+ coorg_send(2,(ispec-1)*5+5) = z2
+ end if
- write(24,*) 'CO'
-
material = kmato(ispec)
mu = elastcoef(2,material)
@@ -470,33 +571,83 @@
! display bad elements that are above 80% of the threshold
if(courant_stability_number >= 0.80 * courant_stability_number_max) then
- write(24,*) '1 0 0 RG GF 0 setgray ST'
+ if ( myrank == 0 ) then
+ write(24,*) '1 0 0 RG GF 0 setgray ST'
+ else
+ RGB_send(ispec) = 1
+ end if
else
! do not color the elements if below the threshold
- write(24,*) 'ST'
+ if ( myrank == 0 ) then
+ write(24,*) 'ST'
+ else
+ RGB_send(ispec) = 0
+ end if
endif
enddo ! end of loop on all the spectral elements
- write(24,*) '%'
- write(24,*) 'grestore'
- write(24,*) 'showpage'
+#ifdef USE_MPI
+ if (myrank == 0 ) then
+
+ do iproc = 1, nproc-1
+ call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ allocate(coorg_recv(2,nspec_recv*5))
+ allocate(RGB_recv(nspec_recv))
+ call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ call MPI_RECV (RGB_recv(1), nspec_recv, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+
+ do ispec = 1, nspec_recv
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',num_ispec
+ write(24,*) 'mark'
+ write(24,681) coorg_recv(1,(ispec-1)*5+1), coorg_recv(2,(ispec-1)*5+1)
+ write(24,681) coorg_recv(1,(ispec-1)*5+2), coorg_recv(2,(ispec-1)*5+2)
+ write(24,681) coorg_recv(1,(ispec-1)*5+3), coorg_recv(2,(ispec-1)*5+3)
+ write(24,681) coorg_recv(1,(ispec-1)*5+4), coorg_recv(2,(ispec-1)*5+4)
+ write(24,681) coorg_recv(1,(ispec-1)*5+5), coorg_recv(2,(ispec-1)*5+5)
+ write(24,*) 'CO'
+ if ( RGB_recv(ispec) == 1) then
+ write(24,*) '1 0 0 RG GF 0 setgray ST'
+ else
+ write(24,*) 'ST'
+ end if
+ end do
+ deallocate(coorg_recv)
+ deallocate(RGB_recv)
+
+ end do
+
+ else
+ call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (RGB_send, nspec, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
- close(24)
+ end if
+#endif
- print *,'End of creation of PostScript file with stability condition'
+ if ( myrank == 0 ) then
+ write(24,*) '%'
+ write(24,*) 'grestore'
+ write(24,*) 'showpage'
+ close(24)
+
+ print *,'End of creation of PostScript file with stability condition'
+ end if
+
!
!--------------------------------------------------------------------------------
!
-
+
+if ( myrank == 0 ) then
print *
print *,'Creating PostScript file with mesh dispersion'
!
!---- open PostScript file
!
- if(any_elastic) then
+ if(any_elastic_glob) then
open(unit=24,file='OUTPUT_FILES/mesh_S_wave_dispersion.ps',status='unknown')
else
open(unit=24,file='OUTPUT_FILES/mesh_P_wave_dispersion.ps',status='unknown')
@@ -575,7 +726,7 @@
write(24,*) '24.35 CM 18.9 CM MV'
write(24,*) usoffset,' CM 2 div neg 0 MR'
write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
- if(any_elastic) then
+ if(any_elastic_glob) then
write(24,*) '(Mesh elastic S-wave dispersion \(red = good, blue = bad\)) show'
else
write(24,*) '(Mesh acoustic P-wave dispersion \(red = good, blue = bad\)) show'
@@ -603,11 +754,16 @@
write(24,*) '% spectral element mesh'
write(24,*) '%'
write(24,*) '0 setgray'
+
+ num_ispec = 0
+ end if
do ispec = 1, nspec
+ if ( myrank == 0 ) then
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',ispec
+ end if
- write(24,*) '% elem ',ispec
-
do i=1,pointsdisp
do j=1,pointsdisp
xinterp(i,j) = 0.d0
@@ -626,8 +782,13 @@
z1 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x1 = x1 * centim
z1 = z1 * centim
- write(24,*) 'mark'
- write(24,681) x1,z1
+ if ( myrank == 0 ) then
+ write(24,*) 'mark'
+ write(24,681) x1,z1
+ else
+ coorg_send(1,(ispec-1)*5+1) = x1
+ coorg_send(2,(ispec-1)*5+1) = z1
+ end if
! draw straight lines if elements have 4 nodes
@@ -636,7 +797,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+2) = x2
+ coorg_send(2,(ispec-1)*5+2) = z2
+ end if
ir=pointsdisp
is=pointsdisp
@@ -644,7 +810,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+3) = x2
+ coorg_send(2,(ispec-1)*5+3) = z2
+ end if
is=pointsdisp
ir=1
@@ -652,7 +823,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+4) = x2
+ coorg_send(2,(ispec-1)*5+4) = z2
+ end if
ir=1
is=2
@@ -660,9 +836,15 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ write(24,*) 'CO'
+ else
+ coorg_send(1,(ispec-1)*5+5) = x2
+ coorg_send(2,(ispec-1)*5+5) = z2
+ end if
- write(24,*) 'CO'
+
material = kmato(ispec)
@@ -717,7 +899,7 @@
distance_max = max(distance_max,distance_max_local)
! display mesh dispersion for S waves if there is at least one elastic element in the mesh
- if(any_elastic) then
+ if(any_elastic_glob) then
! ignore fluid regions with Vs = 0
if(csloc > 0.0001d0) then
@@ -726,20 +908,36 @@
! display very good elements that are above 80% of the threshold in red
if(lambdaS_local >= 0.80 * lambdaSmax) then
- write(24,*) '1 0 0 RG GF 0 setgray ST'
+ if ( myrank == 0 ) then
+ write(24,*) '1 0 0 RG GF 0 setgray ST'
+ else
+ RGB_send(ispec) = 1
+ end if
! display bad elements that are below 120% of the threshold in blue
else if(lambdaS_local <= 1.20 * lambdaSmin) then
- write(24,*) '0 0 1 RG GF 0 setgray ST'
+ if ( myrank == 0 ) then
+ write(24,*) '0 0 1 RG GF 0 setgray ST'
+ else
+ RGB_send(ispec) = 3
+ end if
else
! do not color the elements if not close to the threshold
- write(24,*) 'ST'
+ if ( myrank == 0 ) then
+ write(24,*) 'ST'
+ else
+ RGB_send(ispec) = 0
+ end if
endif
else
! do not color the elements if S-wave velocity undefined
- write(24,*) 'ST'
+ if ( myrank == 0 ) then
+ write(24,*) 'ST'
+ else
+ RGB_send(ispec) = 0
+ end if
endif
! display mesh dispersion for P waves if there is no elastic element in the mesh
@@ -749,33 +947,93 @@
! display very good elements that are above 80% of the threshold in red
if(lambdaP_local >= 0.80 * lambdaPmax) then
- write(24,*) '1 0 0 RG GF 0 setgray ST'
+ if ( myrank == 0 ) then
+ write(24,*) '1 0 0 RG GF 0 setgray ST'
+ else
+ RGB_send(ispec) = 1
+ end if
! display bad elements that are below 120% of the threshold in blue
else if(lambdaP_local <= 1.20 * lambdaPmin) then
- write(24,*) '0 0 1 RG GF 0 setgray ST'
+ if ( myrank == 0 ) then
+ write(24,*) '0 0 1 RG GF 0 setgray ST'
+ else
+ RGB_send(ispec) = 3
+ end if
else
! do not color the elements if not close to the threshold
- write(24,*) 'ST'
+ if ( myrank == 0 ) then
+ write(24,*) 'ST'
+ else
+ RGB_send(ispec) = 0
+ end if
endif
endif
enddo ! end of loop on all the spectral elements
- write(24,*) '%'
- write(24,*) 'grestore'
- write(24,*) 'showpage'
+#ifdef USE_MPI
+ if (myrank == 0 ) then
+
+ do iproc = 1, nproc-1
+ call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ allocate(coorg_recv(2,nspec_recv*5))
+ allocate(RGB_recv(nspec_recv))
+ call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ call MPI_RECV (RGB_recv(1), nspec_recv, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+
+ do ispec = 1, nspec_recv
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',num_ispec
+ write(24,*) 'mark'
+ write(24,681) coorg_recv(1,(ispec-1)*5+1), coorg_recv(2,(ispec-1)*5+1)
+ write(24,681) coorg_recv(1,(ispec-1)*5+2), coorg_recv(2,(ispec-1)*5+2)
+ write(24,681) coorg_recv(1,(ispec-1)*5+3), coorg_recv(2,(ispec-1)*5+3)
+ write(24,681) coorg_recv(1,(ispec-1)*5+4), coorg_recv(2,(ispec-1)*5+4)
+ write(24,681) coorg_recv(1,(ispec-1)*5+5), coorg_recv(2,(ispec-1)*5+5)
+ write(24,*) 'CO'
+ if ( RGB_recv(ispec) == 1) then
+ write(24,*) '1 0 0 RG GF 0 setgray ST'
+ end if
+ if ( RGB_recv(ispec) == 3) then
+ write(24,*) '0 0 1 RG GF 0 setgray ST'
+ end if
+ if ( RGB_recv(ispec) == 0) then
+ write(24,*) 'ST'
+ end if
+
+ end do
+ deallocate(coorg_recv)
+ deallocate(RGB_recv)
- close(24)
+ end do
+
+ else
+ call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (RGB_send, nspec, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
- print *,'End of creation of PostScript file with mesh dispersion'
+ end if
+#endif
+
+ if ( myrank == 0 ) then
+ write(24,*) '%'
+ write(24,*) 'grestore'
+ write(24,*) 'showpage'
+
+ close(24)
+
+ print *,'End of creation of PostScript file with mesh dispersion'
+ end if
+
!
!--------------------------------------------------------------------------------
!
+ if ( myrank == 0 ) then
print *
print *,'Creating PostScript file with velocity model'
@@ -882,10 +1140,14 @@
write(24,*) '%'
write(24,*) '0 setgray'
+ num_ispec = 0
+end if
+
do ispec = 1, nspec
-
- write(24,*) '% elem ',ispec
-
+ if ( myrank == 0 ) then
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',ispec
+ end if
do i=1,pointsdisp
do j=1,pointsdisp
xinterp(i,j) = 0.d0
@@ -904,8 +1166,13 @@
z1 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x1 = x1 * centim
z1 = z1 * centim
- write(24,*) 'mark'
- write(24,681) x1,z1
+ if ( myrank == 0 ) then
+ write(24,*) 'mark'
+ write(24,681) x1,z1
+ else
+ coorg_send(1,(ispec-1)*5+1) = x1
+ coorg_send(2,(ispec-1)*5+1) = z1
+ end if
! draw straight lines if elements have 4 nodes
@@ -914,15 +1181,25 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
-
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+2) = x2
+ coorg_send(2,(ispec-1)*5+2) = z2
+ end if
+
ir=pointsdisp
is=pointsdisp
x2 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+3) = x2
+ coorg_send(2,(ispec-1)*5+3) = z2
+ end if
is=pointsdisp
ir=1
@@ -930,7 +1207,12 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+4) = x2
+ coorg_send(2,(ispec-1)*5+4) = z2
+ end if
ir=1
is=2
@@ -938,10 +1220,14 @@
z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
- write(24,681) x2,z2
-
- write(24,*) 'CO'
-
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ write(24,*) 'CO'
+ else
+ coorg_send(1,(ispec-1)*5+5) = x2
+ coorg_send(2,(ispec-1)*5+5) = z2
+ end if
+
if((vpmax-vpmin)/vpmin > 0.02d0) then
if(assign_external_model) then
! use lower-left corner
@@ -966,18 +1252,331 @@
x1 = 1.d0 - x1
! display P-velocity model using gray levels
- write(24,*) sngl(x1),' setgray GF 0 setgray ST'
-
+ if ( myrank == 0 ) then
+ write(24,*) sngl(x1),' setgray GF 0 setgray ST'
+ else
+ greyscale_send(ispec) = sngl(x1)
+ end if
enddo ! end of loop on all the spectral elements
+#ifdef USE_MPI
+ if (myrank == 0 ) then
+
+ do iproc = 1, nproc-1
+ call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ allocate(coorg_recv(2,nspec_recv*5))
+ allocate(greyscale_recv(nspec_recv))
+ call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ call MPI_RECV (greyscale_recv(1), nspec_recv, MPI_REAL, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+
+ do ispec = 1, nspec_recv
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',num_ispec
+ write(24,*) 'mark'
+ write(24,681) coorg_recv(1,(ispec-1)*5+1), coorg_recv(2,(ispec-1)*5+1)
+ write(24,681) coorg_recv(1,(ispec-1)*5+2), coorg_recv(2,(ispec-1)*5+2)
+ write(24,681) coorg_recv(1,(ispec-1)*5+3), coorg_recv(2,(ispec-1)*5+3)
+ write(24,681) coorg_recv(1,(ispec-1)*5+4), coorg_recv(2,(ispec-1)*5+4)
+ write(24,681) coorg_recv(1,(ispec-1)*5+5), coorg_recv(2,(ispec-1)*5+5)
+ write(24,*) 'CO'
+ write(24,*) greyscale_recv(ispec), ' setgray GF 0 setgray ST'
+
+ end do
+ deallocate(coorg_recv)
+ deallocate(greyscale_recv)
+
+ end do
+
+ else
+ call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (greyscale_send, nspec, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
+
+ end if
+#endif
+
+ if ( myrank == 0 ) then
+ write(24,*) '%'
+ write(24,*) 'grestore'
+ write(24,*) 'showpage'
+
+ close(24)
+
+ print *,'End of creation of PostScript file with velocity model'
+
+ end if
+
+
+ print *
+ print *,'Creating PostScript file with partitioning'
+
+ if ( myrank == 0 ) then
+!
+!---- open PostScript file
+!
+ open(unit=24,file='OUTPUT_FILES/mesh_partition.ps',status='unknown')
+
+!
+!---- write PostScript header
+!
+ write(24,10) simulation_title
+ write(24,*) '/CM {28.5 mul} def'
+ write(24,*) '/LR {rlineto} def'
+ write(24,*) '/LT {lineto} def'
+ write(24,*) '/L {lineto} def'
+ write(24,*) '/MR {rmoveto} def'
+ write(24,*) '/MV {moveto} def'
+ write(24,*) '/M {moveto} def'
+ write(24,*) '/ST {stroke} def'
+ write(24,*) '/CP {closepath} def'
+ write(24,*) '/RG {setrgbcolor} def'
+ write(24,*) '/GF {gsave fill grestore} def'
+ write(24,*) '% different useful symbols'
+ write(24,*) '/Point {2 0 360 arc CP 0 setgray fill} def'
+ write(24,*) '/VDot {-0.75 -1.5 MR 1.5 0 LR 0 3. LR -1.5 0 LR'
+ write(24,*) 'CP fill} def'
+ write(24,*) '/HDot {-1.5 -0.75 MR 3. 0 LR 0 1.5 LR -3. 0 LR'
+ write(24,*) 'CP fill} def'
+ write(24,*) '/Cross {gsave 0.05 CM setlinewidth'
+ write(24,*) 'gsave 3 3 MR -6. -6. LR ST grestore'
+ write(24,*) 'gsave 3 -3 MR -6. 6. LR ST grestore'
+ write(24,*) '0.01 CM setlinewidth} def'
+ write(24,*) '/SmallLine {MV 0.07 CM 0 rlineto} def'
+ write(24,*) '/Diamond {gsave 0.05 CM setlinewidth 0 4.2 MR'
+ write(24,*) '-3 -4.2 LR 3 -4.2 LR 3 4.2 LR CP ST'
+ write(24,*) 'grestore 0.01 CM setlinewidth} def'
write(24,*) '%'
+ write(24,*) '% macro to draw the contour of the elements'
+ write(24,*) '/CO {M counttomark 2 idiv {L} repeat cleartomark CP} def'
+ write(24,*) '%'
+ write(24,*) '.01 CM setlinewidth'
+ write(24,*) '/Times-Roman findfont'
+ write(24,*) '.35 CM scalefont setfont'
+ write(24,*) '%'
+ write(24,*) '/vshift ',-height/2,' CM def'
+ write(24,*) '/Rshow { currentpoint stroke MV'
+ write(24,*) 'dup stringwidth pop neg vshift MR show } def'
+ write(24,*) '/Cshow { currentpoint stroke MV'
+ write(24,*) 'dup stringwidth pop -2 div vshift MR show } def'
+ write(24,*) '/fN {/Helvetica-Bold findfont ',height,' CM scalefont setfont} def'
+ write(24,*) '%'
+ write(24,*) 'gsave newpath 90 rotate'
+ write(24,*) '0 ',-sizez,' CM translate 1. 1. scale'
+ write(24,*) '%'
+
+!
+!--- write captions of PostScript figure
+!
+ write(24,*) '0 setgray'
+ write(24,*) '/Times-Roman findfont'
+ write(24,*) '.5 CM scalefont setfont'
+
+ write(24,*) '%'
+ write(24,*) '/Times-Roman findfont'
+ write(24,*) '.6 CM scalefont setfont'
+ write(24,*) '.4 .9 .9 setrgbcolor'
+ write(24,*) '11 CM 1.1 CM MV'
+ write(24,*) '(X axis) show'
+ write(24,*) '%'
+ write(24,*) '1.4 CM 9.5 CM MV'
+ write(24,*) 'currentpoint gsave translate 90 rotate 0 0 moveto'
+ write(24,*) '(Y axis) show'
write(24,*) 'grestore'
- write(24,*) 'showpage'
+ write(24,*) '%'
+ write(24,*) '/Times-Roman findfont'
+ write(24,*) '.7 CM scalefont setfont'
+ write(24,*) '.8 0 .8 setrgbcolor'
+ write(24,*) '24.35 CM 18.9 CM MV'
+ write(24,*) usoffset,' CM 2 div neg 0 MR'
+ write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
+ write(24,*) '(Mesh stability condition \(red = bad\)) show'
+ write(24,*) 'grestore'
+ write(24,*) '25.35 CM 18.9 CM MV'
+ write(24,*) usoffset,' CM 2 div neg 0 MR'
+ write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
+ write(24,*) '(',simulation_title,') show'
+ write(24,*) 'grestore'
+ write(24,*) '26.45 CM 18.9 CM MV'
+ write(24,*) usoffset,' CM 2 div neg 0 MR'
+ write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
+ write(24,*) '(2D Spectral Element Method) show'
+ write(24,*) 'grestore'
- close(24)
+ write(24,*) '%'
+ write(24,*) '1 1 scale'
+ write(24,*) '%'
- print *,'End of creation of PostScript file with velocity model'
+!
+!---- draw the spectral element mesh
+!
+ write(24,*) '%'
+ write(24,*) '% spectral element mesh'
+ write(24,*) '%'
+ write(24,*) '0 setgray'
+ num_ispec = 0
+ end if
+
+ do ispec = 1, nspec
+
+ if ( myrank == 0 ) then
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',ispec
+ end if
+
+ do i=1,pointsdisp
+ do j=1,pointsdisp
+ xinterp(i,j) = 0.d0
+ zinterp(i,j) = 0.d0
+ do in = 1,ngnod
+ nnum = knods(in,ispec)
+ xinterp(i,j) = xinterp(i,j) + shapeint(in,i,j)*coorg(1,nnum)
+ zinterp(i,j) = zinterp(i,j) + shapeint(in,i,j)*coorg(2,nnum)
+ enddo
+ enddo
+ enddo
+
+ is = 1
+ ir = 1
+ x1 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
+ z1 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
+ x1 = x1 * centim
+ z1 = z1 * centim
+ if ( myrank == 0 ) then
+ write(24,*) 'mark'
+ write(24,681) x1,z1
+ else
+ coorg_send(1,(ispec-1)*5+1) = x1
+ coorg_send(2,(ispec-1)*5+1) = z1
+ end if
+
+! draw straight lines if elements have 4 nodes
+
+ ir=pointsdisp
+ x2 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+2) = x2
+ coorg_send(2,(ispec-1)*5+2) = z2
+ end if
+
+ ir=pointsdisp
+ is=pointsdisp
+ x2 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+3) = x2
+ coorg_send(2,(ispec-1)*5+3) = z2
+ end if
+
+ is=pointsdisp
+ ir=1
+ x2 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ else
+ coorg_send(1,(ispec-1)*5+4) = x2
+ coorg_send(2,(ispec-1)*5+4) = z2
+ end if
+
+ ir=1
+ is=2
+ x2 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*ratio_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ if ( myrank == 0 ) then
+ write(24,681) x2,z2
+ write(24,*) 'CO'
+ else
+ coorg_send(1,(ispec-1)*5+5) = x2
+ coorg_send(2,(ispec-1)*5+5) = z2
+ end if
+
+
+ if ( myrank == 0 ) then
+ write(24,*) '0.3 0.3 0.4 RG GF 0 setgray ST'
+ end if
+
+ enddo ! end of loop on all the spectral elements
+
+#ifdef USE_MPI
+ if (myrank == 0 ) then
+
+ RGB_step = 0.6 / real(nproc_per_cat)
+ nproc_per_cat = ceiling(real(nproc)/real(3))
+
+ do iproc = 1, nproc-1
+
+ color_cat = ceiling(real(iproc)/real(nproc_per_cat))
+ if ( color_cat == 1 ) then
+ R_val = 0.2 + RGB_step * modulo(iproc,nproc_per_cat)
+ G_val = 1. - R_val
+ B_val = 0
+ end if
+ if ( color_cat == 2 ) then
+ R_val = 0
+ G_val = 0.2 + RGB_step * modulo(iproc,nproc_per_cat)
+ B_val = 1. - G_val
+ end if
+ if ( color_cat == 3 ) then
+ R_val = 0.2 + RGB_step * modulo(iproc,nproc_per_cat)
+ G_val = 0
+ B_val = 1. - R_val
+ end if
+
+ call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ allocate(coorg_recv(2,nspec_recv*5))
+ call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+
+ do ispec = 1, nspec_recv
+ num_ispec = num_ispec + 1
+ write(24,*) '% elem ',num_ispec
+ write(24,*) 'mark'
+ write(24,681) coorg_recv(1,(ispec-1)*5+1), coorg_recv(2,(ispec-1)*5+1)
+ write(24,681) coorg_recv(1,(ispec-1)*5+2), coorg_recv(2,(ispec-1)*5+2)
+ write(24,681) coorg_recv(1,(ispec-1)*5+3), coorg_recv(2,(ispec-1)*5+3)
+ write(24,681) coorg_recv(1,(ispec-1)*5+4), coorg_recv(2,(ispec-1)*5+4)
+ write(24,681) coorg_recv(1,(ispec-1)*5+5), coorg_recv(2,(ispec-1)*5+5)
+ write(24,*) 'CO'
+
+ write(24,*) R_val, G_val, B_val, ' RG GF 0 setgray ST'
+
+ end do
+ deallocate(coorg_recv)
+
+ end do
+
+ else
+ call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
+ call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
+
+ end if
+#endif
+
+ if ( myrank == 0 ) then
+ write(24,*) '%'
+ write(24,*) 'grestore'
+ write(24,*) 'showpage'
+
+ close(24)
+
+ print *,'End of creation of PostScript file with partitioning'
+ end if
+
+
+
10 format('%!PS-Adobe-2.0',/,'%%',/,'%% Title: ',a50,/,'%% Created by: Specfem2D',/,'%% Author: Dimitri Komatitsch',/,'%%')
681 format(f6.2,1x,f6.2)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -12,7 +12,7 @@
!========================================================================
subroutine compute_forces_acoustic(npoin,nspec,nelemabs,numat, &
- iglob_source,ispec_selected_source,source_type,it,NSTEP,anyabs, &
+ iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs, &
assign_external_model,initialfield,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
@@ -27,7 +27,7 @@
include "constants.h"
- integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,source_type,it,NSTEP
+ integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP
logical :: anyabs,assign_external_model,initialfield
@@ -334,22 +334,23 @@
! --- add the source
if(.not. initialfield) then
+ if (is_proc_source == 1 ) then
! collocated force
! beware, for acoustic medium, source is a pressure source
- if(source_type == 1) then
- if(.not. elastic(ispec_selected_source)) then
- potential_dot_dot_acoustic(iglob_source) = potential_dot_dot_acoustic(iglob_source) + source_time_function(it)
- endif
+ if(source_type == 1) then
+ if(.not. elastic(ispec_selected_source)) then
+ potential_dot_dot_acoustic(iglob_source) = potential_dot_dot_acoustic(iglob_source) + source_time_function(it)
+ endif
! moment tensor
- else if(source_type == 2) then
+ else if(source_type == 2) then
- if(.not. elastic(ispec_selected_source)) stop 'cannot have moment tensor source in acoustic element'
+ if(.not. elastic(ispec_selected_source)) stop 'cannot have moment tensor source in acoustic element'
- endif
-
+ endif
+ end if
else
- stop 'wrong source type'
+ stop 'wrong source type'
endif
end subroutine compute_forces_acoustic
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -12,7 +12,7 @@
!========================================================================
subroutine compute_forces_elastic(npoin,nspec,nelemabs,numat,iglob_source, &
- ispec_selected_source,source_type,it,NSTEP,anyabs,assign_external_model, &
+ ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic,density,elastcoef,xix,xiz,gammax,gammaz, &
@@ -27,7 +27,7 @@
include "constants.h"
- integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,source_type,it,NSTEP
+ integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP
logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
@@ -443,29 +443,31 @@
! --- add the source
if(.not. initialfield) then
+ if ( is_proc_source == 1 ) then
! collocated force
! beware, for acoustic medium, source is a potential, therefore source time function
! gives shape of velocity, not displacement
- if(source_type == 1) then
- if(elastic(ispec_selected_source)) then
- accel_elastic(1,iglob_source) = accel_elastic(1,iglob_source) - sin(angleforce)*source_time_function(it)
- accel_elastic(2,iglob_source) = accel_elastic(2,iglob_source) + cos(angleforce)*source_time_function(it)
- endif
+ if(source_type == 1) then
+ if(elastic(ispec_selected_source)) then
+ accel_elastic(1,iglob_source) = accel_elastic(1,iglob_source) - sin(angleforce)*source_time_function(it)
+ accel_elastic(2,iglob_source) = accel_elastic(2,iglob_source) + cos(angleforce)*source_time_function(it)
+ endif
! moment tensor
- else if(source_type == 2) then
+ else if(source_type == 2) then
- if(elastic(ispec_selected_source)) then
+ if(elastic(ispec_selected_source)) then
! add source array
- do j=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,ispec_selected_source)
- accel_elastic(:,iglob) = accel_elastic(:,iglob) + sourcearray(:,i,j)*source_time_function(it)
- enddo
- enddo
- endif
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source)
+ accel_elastic(:,iglob) = accel_elastic(:,iglob) + sourcearray(:,i,j)*source_time_function(it)
+ enddo
+ enddo
+ endif
- endif
+ endif
+ end if
else
stop 'wrong source type'
Added: seismo/2D/SPECFEM2D/trunk/constants_unstruct.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants_unstruct.h 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/constants_unstruct.h 2007-12-07 23:53:26 UTC (rev 8524)
@@ -0,0 +1,9 @@
+
+! Number of nodes per elements.
+integer, parameter :: ESIZE = 4
+
+! Max number of neighbours per elements.
+ integer :: max_neighbour=30
+
+! Max number of elements that can contain the same node.
+ integer :: nsize=20
Modified: seismo/2D/SPECFEM2D/trunk/create_color_image.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -224,8 +224,8 @@
close(27)
! open image file and create system command to convert image to more convenient format
- write(system_command,"('cd OUTPUT_FILES ; convert image',i6.6,'.pnm image',i6.6,'.gif ; rm image',i6.6,'.pnm')") it,it,it
-
+ write(system_command,"('cd OUTPUT_FILES ; convert image',i6.6,'.pnm image',i6.6,'.gif')") it,it
+
! call the system to convert image to GIF
call system(system_command)
Modified: seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -12,7 +12,7 @@
!========================================================================
subroutine enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,ispecnum_acoustic_surface,iedgenum_acoustic_surface, &
+ potential_acoustic,acoustic_surface, &
ibool,nelem_acoustic_surface,npoin,nspec)
! free surface for an acoustic medium
@@ -25,7 +25,7 @@
integer :: nelem_acoustic_surface,npoin,nspec
- integer, dimension(nelem_acoustic_surface) :: ispecnum_acoustic_surface,iedgenum_acoustic_surface
+ integer, dimension(5,nelem_acoustic_surface) :: acoustic_surface
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
@@ -36,39 +36,21 @@
!---
integer :: ispec_acoustic_surface,ispec,iedge,i,j,iglob
-
- do ispec_acoustic_surface = 1,nelem_acoustic_surface
-
- ispec = ispecnum_acoustic_surface(ispec_acoustic_surface)
- iedge = iedgenum_acoustic_surface(ispec_acoustic_surface)
-
- if(iedge == IBOTTOM .or. iedge == ITOP) then
- if(iedge == IBOTTOM) then
- j = 1
- else
- j = NGLLZ
- endif
- do i = 1,NGLLX
- iglob = ibool(i,j,ispec)
- potential_acoustic(iglob) = ZERO
- potential_dot_acoustic(iglob) = ZERO
- potential_dot_dot_acoustic(iglob) = ZERO
- enddo
- else
- if(iedge == ILEFT) then
- i = 1
- else
- i = NGLLX
- endif
- do j = 1,NGLLZ
- iglob = ibool(i,j,ispec)
- potential_acoustic(iglob) = ZERO
- potential_dot_acoustic(iglob) = ZERO
- potential_dot_dot_acoustic(iglob) = ZERO
- enddo
- endif
-
- enddo
-
+
+ do ispec_acoustic_surface = 1, nelem_acoustic_surface
+
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+ do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
+ do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
+ iglob = ibool(i,j,ispec)
+ potential_acoustic(iglob) = ZERO
+ potential_dot_acoustic(iglob) = ZERO
+ potential_dot_dot_acoustic(iglob) = ZERO
+
+ end do
+ end do
+
+ end do
+
end subroutine enforce_acoustic_free_surface
Modified: seismo/2D/SPECFEM2D/trunk/locate_receivers.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_receivers.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/locate_receivers.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -15,14 +15,19 @@
!---- locate_receivers finds the correct position of the receivers
!----
- subroutine locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,st_xval,st_zval,ispec_selected_rec, &
- xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo)
+ subroutine locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank, &
+ st_xval,st_zval,ispec_selected_rec, &
+ xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo)
implicit none
include "constants.h"
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
integer nrec,nspec,npoin,ngnod,npgeo
+ integer, intent(in) :: nproc, myrank
integer knods(ngnod,nspec)
double precision coorg(NDIM,npgeo)
@@ -34,7 +39,8 @@
integer nrec_dummy,irec,i,j,ispec,iglob,iter_loop,ix_initial_guess,iz_initial_guess
- double precision x_source,z_source,dist,stele,stbur,distance_receiver
+ double precision x_source,z_source,dist,stele,stbur
+ double precision, dimension(nrec) :: distance_receiver
double precision xi,gamma,dx,dz,dxi,dgamma
! Gauss-Lobatto-Legendre points of integration
@@ -44,11 +50,12 @@
double precision x,z,xix,xiz,gammax,gammaz,jacobian
! use dynamic allocation
- double precision distmin
+ double precision distmin, dist_glob
double precision, dimension(:), allocatable :: final_distance
! receiver information
- integer, dimension(nrec) :: ispec_selected_rec
+ integer :: nrecloc, is_proc_receiver, nb_proc_receiver
+ integer, dimension(nrec) :: ispec_selected_rec, recloc
double precision, dimension(nrec) :: xi_receiver,gamma_receiver
! station information for writing the seismograms
@@ -57,6 +64,14 @@
double precision, dimension(nrec) :: st_xval,st_zval
+
+ double precision, dimension(nrec,nproc) :: gather_final_distance
+ double precision, dimension(nrec,nproc) :: gather_xi_receiver, gather_gamma_receiver
+ integer, dimension(nrec,nproc) :: gather_ispec_selected_rec
+ integer, dimension(nrec), intent(inout) :: which_proc_receiver
+ integer :: ierror
+
+
! **************
write(IOUT,*)
@@ -68,7 +83,7 @@
write(IOUT,*)
! get number of stations from receiver file
- open(unit=1,file='DATA/STATIONS',status='old')
+ open(unit=1,file='DATA/STATIONS',status='old',action='read')
read(1,*) nrec_dummy
if(nrec_dummy /= nrec) stop 'problem with number of receivers'
@@ -88,7 +103,7 @@
if(abs(stbur) > TINYVAL) stop 'stations with non-zero burial not implemented yet'
! compute distance between source and receiver
- distance_receiver = sqrt((st_zval(irec)-z_source)**2 + (st_xval(irec)-x_source)**2)
+ distance_receiver(irec) = sqrt((st_zval(irec)-z_source)**2 + (st_xval(irec)-x_source)**2)
do ispec=1,nspec
@@ -114,6 +129,7 @@
! end of loop on all the spectral elements
enddo
+
! ****************************************
! find the best (xi,gamma) for each receiver
! ****************************************
@@ -163,33 +179,93 @@
! compute final distance between asked and found
final_distance(irec) = sqrt((st_xval(irec)-x)**2 + (st_zval(irec)-z)**2)
+enddo
+
+! close receiver file
+close(1)
+
+! elect one process for each receiver.
+#ifdef USE_MPI
+call MPI_GATHER(final_distance(1),nrec,MPI_DOUBLE_PRECISION,&
+ gather_final_distance(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
+call MPI_GATHER(xi_receiver(1),nrec,MPI_DOUBLE_PRECISION,&
+ gather_xi_receiver(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
+call MPI_GATHER(gamma_receiver(1),nrec,MPI_DOUBLE_PRECISION,&
+ gather_gamma_receiver(1,1),nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierror)
+call MPI_GATHER(ispec_selected_rec(1),nrec,MPI_INTEGER,&
+ gather_ispec_selected_rec(1,1),nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
+
+if ( myrank == 0 ) then
+ do irec = 1, nrec
+ which_proc_receiver(irec:irec) = minloc(gather_final_distance(irec,:)) - 1
+
+ end do
+end if
+
+call MPI_BCAST(which_proc_receiver(1),nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ierror)
+
+#else
+print *, 'POY', nproc
+call flush(6)
+gather_final_distance(:,1) = final_distance(:)
+print *, 'POY'
+call flush(6)
+gather_xi_receiver(:,1) = xi_receiver(:)
+gather_gamma_receiver(:,1) = gamma_receiver(:)
+gather_ispec_selected_rec(:,1) = ispec_selected_rec(:)
+
+which_proc_receiver(:) = 0
+
+#endif
+
+nrecloc = 0
+do irec = 1, nrec
+ if ( which_proc_receiver(irec) == myrank ) then
+ nrecloc = nrecloc + 1
+ recloc(nrecloc) = irec
+ end if
+
+end do
+
+
+if ( myrank == 0 ) then
+
+ do irec = 1, nrec
write(IOUT,*)
write(IOUT,*) 'Station # ',irec,' ',station_name(irec),network_name(irec)
- if(final_distance(irec) == HUGEVAL) stop 'error locating receiver'
+ if(gather_final_distance(irec,which_proc_receiver(irec)+1) == HUGEVAL) stop 'error locating receiver'
write(IOUT,*) ' original x: ',sngl(st_xval(irec))
write(IOUT,*) ' original z: ',sngl(st_zval(irec))
- write(IOUT,*) ' distance from source: ',sngl(distance_receiver)
- write(IOUT,*) 'closest estimate found: ',sngl(final_distance(irec)),' m away'
- write(IOUT,*) ' in element ',ispec_selected_rec(irec)
- write(IOUT,*) ' at xi,gamma coordinates = ',xi_receiver(irec),gamma_receiver(irec)
+ write(IOUT,*) ' distance from source: ',sngl(distance_receiver(irec))
+ write(IOUT,*) 'closest estimate found: ',sngl(gather_final_distance(irec,which_proc_receiver(irec)+1)),' m away'
+ write(IOUT,*) ' in element ',gather_ispec_selected_rec(irec,which_proc_receiver(irec)+1)
+ write(IOUT,*) ' at process ', which_proc_receiver(irec)
+ write(IOUT,*) ' at xi,gamma coordinates = ',gather_xi_receiver(irec,which_proc_receiver(irec)+1),&
+ gather_gamma_receiver(irec,which_proc_receiver(irec)+1)
write(IOUT,*)
- enddo
+ end do
-! close receiver file
- close(1)
! display maximum error for all the receivers
- write(IOUT,*) 'maximum error in location of all the receivers: ',sngl(maxval(final_distance(:))),' m'
+ !write(IOUT,*) 'maximum error in location of all the receivers: ',sngl(maxval(final_distance(:))),' m'
write(IOUT,*)
write(IOUT,*) 'end of receiver detection'
write(IOUT,*)
+
+end if
+
! deallocate arrays
deallocate(final_distance)
+#ifdef USE_MPI
+ call MPI_BARRIER(MPI_COMM_WORLD,ierror)
+#endif
+
+
end subroutine locate_receivers
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -12,7 +12,7 @@
!========================================================================
subroutine locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,source_type,ix_source,iz_source, &
- ispec_source,iglob_source)
+ ispec_source,iglob_source,is_proc_source,nb_proc_source)
!
!----- calculer la position reelle de la source
@@ -21,6 +21,9 @@
implicit none
include "constants.h"
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
integer npoin,nspec,source_type
integer ibool(NGLLX,NGLLZ,nspec)
@@ -28,12 +31,17 @@
double precision x_source,z_source
double precision coord(NDIM,npoin)
+ integer, intent(inout) :: is_proc_source, nb_proc_source
+
integer ip,ix,iz,numelem,ilowx,ilowz,ihighx,ihighz,ix_source,iz_source,ispec_source,iglob_source
- double precision distminmax,distmin,xp,zp,dist
+ double precision distminmax,distmin,xp,zp,dist, dist_glob
- write(iout,200)
+ integer :: ierror
+
+ is_proc_source = 0
+
distminmax = -HUGEVAL
distmin = +HUGEVAL
@@ -42,6 +50,10 @@
ilowz = 1
ihighx = NGLLX
ihighz = NGLLZ
+!!$ ilowx = 2
+!!$ ilowz = 2
+!!$ ihighx = NGLLX-1
+!!$ ihighz = NGLLZ-1
! on ne fait la recherche que sur l'interieur de l'element si source explosive
if(source_type == 2) then
@@ -80,14 +92,54 @@
distminmax = max(distmin,distminmax)
- write(iout,"(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3)") x_source,z_source,coord(1,iglob_source),coord(2,iglob_source),distmin
- write(iout,*)
- write(iout,*)
- write(iout,"('Maximum distance between asked and real =',f12.3)") distminmax
+#ifdef USE_MPI
+ ! global minimum distance computed over all processes
+ call MPI_ALLREDUCE (distminmax, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
+
+#else
+ dist_glob = distminmax
+
+#endif
+
+
+ ! check if this process contains the source
+ if ( dist_glob == distminmax ) then
+ is_proc_source = 1
+ end if
+
+
+#ifdef USE_MPI
+ ! determining the number of processes that contain the source (useful when the source is located on an interface)
+ call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
+
+#else
+ nb_proc_source = is_proc_source
+
+#endif
+
+ if ( nb_proc_source < 1 ) then
+ stop "error locating force source"
+ end if
+
+ if ( is_proc_source == 1 ) then
+ write(iout,200)
+
+ write(iout,"(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3,1x,i5.5)") x_source,z_source, &
+ coord(1,iglob_source),coord(2,iglob_source),distmin,nb_proc_source
+ write(iout,*)
+ write(iout,*)
+ write(iout,"('Maximum distance between asked and real =',f12.3)") distminmax
+
+ end if
+
+#ifdef USE_MPI
+ call MPI_BARRIER(MPI_COMM_WORLD,ierror)
+#endif
+
200 format(//1x,48('=')/,' = S o u r c e s ', &
'r e a l p o s i t i o n s ='/1x,48('=')// &
- ' Source x-asked z-asked x-obtain z-obtain dist'/)
+ ' Source x-asked z-asked x-obtain z-obtain dist nb_proc_source'/)
end subroutine locate_source_force
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -16,11 +16,14 @@
!----
subroutine locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
- ispec_selected_source,xi_source,gamma_source,coorg,knods,ngnod,npgeo)
+ ispec_selected_source,is_proc_source,nb_proc_source,nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod,npgeo)
implicit none
include "constants.h"
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
integer nspec,npoin,ngnod,npgeo
@@ -42,12 +45,17 @@
double precision zigll(NGLLZ)
double precision x,z,xix,xiz,gammax,gammaz,jacobian
- double precision distmin,final_distance
+ double precision distmin,final_distance,dist_glob
! source information
- integer ispec_selected_source
+ integer ispec_selected_source,is_proc_source,nb_proc_source
+ integer, intent(in) :: nproc, myrank
double precision xi_source,gamma_source
+ integer, dimension(1:nproc) :: allgather_is_proc_source
+ integer, dimension(1) :: locate_is_proc_source
+ integer :: ierror
+
! **************
write(IOUT,*)
@@ -59,30 +67,75 @@
! set distance to huge initial value
distmin=HUGEVAL
- do ispec=1,nspec
+ is_proc_source = 0
+ do ispec=1,nspec
+
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
- do j=2,NGLLZ-1
- do i=2,NGLLX-1
+ do j=2,NGLLZ-1
+ do i=2,NGLLX-1
- iglob = ibool(i,j,ispec)
- dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
-
-! keep this point if it is closer to the source
- if(dist < distmin) then
+ iglob = ibool(i,j,ispec)
+ dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
+
+! keep this point if it is closer to the source
+ if(dist < distmin) then
distmin = dist
ispec_selected_source = ispec
ix_initial_guess = i
iz_initial_guess = j
- endif
-
- enddo
+ endif
+
enddo
+ enddo
! end of loop on all the spectral elements
- enddo
+ enddo
+
+#ifdef USE_MPI
+ ! global minimum distance computed over all processes
+ call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
+
+#else
+ dist_glob = distmin
+
+#endif
+
+
+ ! check if this process contains the source
+ if ( dist_glob == distmin ) then
+ is_proc_source = 1
+ end if
+
+
+#ifdef USE_MPI
+ ! determining the number of processes that contain the source (useful when the source is located on an interface)
+ call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
+
+#else
+ nb_proc_source = is_proc_source
+
+#endif
+
+
+#ifdef USE_MPI
+ ! when several processes contain the source, we elect one of them (minimum rank).
+ if ( nb_proc_source > 1 ) then
+
+ call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
+ locate_is_proc_source = maxloc(allgather_is_proc_source) - 1
+
+ if ( myrank /= locate_is_proc_source(1) ) then
+ is_proc_source = 0
+ end if
+ nb_proc_source = 1
+
+ end if
+
+#endif
+
! ****************************************
! find the best (xi,gamma) for each source
! ****************************************
@@ -132,21 +185,27 @@
! compute final distance between asked and found
final_distance = sqrt((x_source-x)**2 + (z_source-z)**2)
- write(IOUT,*)
- write(IOUT,*) 'Moment-tensor source:'
+ if ( is_proc_source == 1 ) then
+ write(IOUT,*)
+ write(IOUT,*) 'Moment-tensor source:'
+
+ if(final_distance == HUGEVAL) stop 'error locating moment-tensor source'
+
+ write(IOUT,*) ' original x: ',sngl(x_source)
+ write(IOUT,*) ' original z: ',sngl(z_source)
+ write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
+ write(IOUT,*) ' in element ',ispec_selected_source
+ write(IOUT,*) ' at xi,gamma coordinates = ',xi_source,gamma_source
+ write(IOUT,*)
- if(final_distance == HUGEVAL) stop 'error locating moment-tensor source'
+ write(IOUT,*)
+ write(IOUT,*) 'end of moment-tensor source detection'
+ write(IOUT,*)
+ end if
- write(IOUT,*) ' original x: ',sngl(x_source)
- write(IOUT,*) ' original z: ',sngl(z_source)
- write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
- write(IOUT,*) ' in element ',ispec_selected_source
- write(IOUT,*) ' at xi,gamma coordinates = ',xi_source,gamma_source
- write(IOUT,*)
+#ifdef USE_MPI
+ call MPI_BARRIER(MPI_COMM_WORLD,ierror)
+#endif
- write(IOUT,*)
- write(IOUT,*) 'end of moment-tensor source detection'
- write(IOUT,*)
-
end subroutine locate_source_moment_tensor
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -38,8 +38,9 @@
! number=2,
! pages={368-392}}
- program meshfem2D
+program meshfem2D
+ use part_unstruct
implicit none
include "constants.h"
@@ -54,7 +55,7 @@
! to store density and velocity model
double precision, dimension(:), allocatable :: rho,cp,cs,aniso3,aniso4
integer, dimension(:), allocatable :: icodemat
- integer, dimension(:,:), allocatable :: num_material
+ integer, dimension(:), allocatable :: num_material
! interface data
integer interface_current,ipoint_current,number_of_interfaces,npoints_interface_bottom,npoints_interface_top
@@ -90,13 +91,13 @@
double precision, dimension(:), allocatable :: xdeb,zdeb,xfin,zfin
logical interpol,gnuplot,assign_external_model,outputgrid
- logical abstop,absbottom,absleft,absright
+ logical abstop,absbottom,absleft,absright,any_abs
logical source_surf,meshvect,initialfield,modelvect,boundvect
logical TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
logical, dimension(:), allocatable :: enreg_surf
- integer, external :: num
+ integer, external :: num_4, num_9
double precision, external :: value_spline
! flag to indicate an anisotropic material
@@ -108,6 +109,67 @@
! ignore variable name field (junk) at the beginning of each input line
logical, parameter :: IGNORE_JUNK = .true.,DONT_IGNORE_JUNK = .false.
+! parameters for external mesh
+ logical :: read_external_mesh
+ character(len=256) :: mesh_file, nodes_coords_file, materials_file, free_surface_file, absorbing_surface_file
+
+! variables used for storing info about the mesh and partitions
+ integer, dimension(:), pointer :: elmnts
+ integer, dimension(:), pointer :: elmnts_bis
+ double precision, dimension(:,:), pointer :: nodes_coords
+ integer, dimension(:,:), pointer :: acoustic_surface
+ integer, dimension(:,:), pointer :: abs_surface
+
+ integer, dimension(:), pointer :: xadj
+ integer, dimension(:), pointer :: adjncy
+ integer, dimension(:), pointer :: nnodes_elmnts
+ integer, dimension(:), pointer :: nodes_elmnts
+
+ integer, dimension(:), pointer :: vwgt
+ integer, dimension(:), pointer :: adjwgt
+ integer, dimension(:), pointer :: part
+
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, dimension(:), pointer :: tab_size_interfaces, tab_interfaces
+ integer, dimension(:), allocatable :: my_interfaces
+ integer, dimension(:), allocatable :: my_nb_interfaces
+
+ integer :: nedges_coupled, nedges_coupled_loc
+ integer, dimension(:,:), pointer :: edges_coupled
+
+ integer :: num_start
+ integer :: nelmnts
+ integer :: num_elmnt
+ integer :: nnodes
+ integer :: num_node
+ integer :: nb_edges
+ integer :: ninterfaces
+ integer :: my_ninterface
+ integer :: nelem_acoustic_surface_loc, nelemabs_loc
+ logical, dimension(:,:), pointer :: abs_surface_char
+ integer, dimension(:), pointer :: abs_surface_merge
+ integer :: nelemabs_merge
+ integer, dimension(:), pointer :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right
+
+! variables used for partitionning
+ integer :: nproc
+ integer :: partitionning_method
+ character(len=256) :: partitionning_strategy
+ character(len=256) :: scotch_strategy
+ integer, dimension(0:4) :: metis_options
+ character(len=256) :: prname
+
+ integer :: edgecut
+ integer :: iproc
+
+ ! variable de test
+ integer :: aaa
+
+
! ***
! *** read the parameter file
! ***
@@ -125,83 +187,184 @@
write(*,*) title
print *
+! read info about external mesh
+ call read_value_logical(IIN,IGNORE_JUNK,read_external_mesh)
+ call read_value_string(IIN,IGNORE_JUNK,mesh_file)
+ call read_value_string(IIN,IGNORE_JUNK,nodes_coords_file)
+ call read_value_string(IIN,IGNORE_JUNK,materials_file)
+ call read_value_string(IIN,IGNORE_JUNK,free_surface_file)
+ call read_value_string(IIN,IGNORE_JUNK,absorbing_surface_file)
+
+! read info about partitionning
+ call read_value_integer(IIN,IGNORE_JUNK,nproc)
+ if ( nproc <= 0 ) then
+ print *, 'Number of processes (nproc) must be greater than or equal to one.'
+ stop
+ end if
+
+#ifndef USE_MPI
+ if ( nproc > 1 ) then
+ print *, 'Number of processes (nproc) must be equal to one when not using MPI.'
+ print *, 'Please recompile with -DUSE_MPI in order to enable use of MPI.'
+ stop
+ end if
+
+#endif
+
+ call read_value_integer(IIN,IGNORE_JUNK,partitionning_method)
+ call read_value_string(IIN,IGNORE_JUNK,partitionning_strategy)
+ select case(partitionning_method)
+ case(1)
+ case(2)
+ partitionning_strategy = trim(partitionning_strategy)
+ if ( partitionning_strategy(1:1) == '0' ) then
+ metis_options(0) = 0
+ else
+ do i = 1, 5
+ metis_options = iachar(partitionning_strategy(i:i)) - iachar('0')
+ end do
+ end if
+
+ case(3)
+ scotch_strategy = trim(partitionning_strategy)
+
+ case default
+ print *, 'Invalid partionning method number.'
+ print *, 'Partionning method', partitionning_method, 'was requested, but is not available.'
+ stop
+ end select
+
+
! read grid parameters
call read_value_double_precision(IIN,IGNORE_JUNK,xmin)
call read_value_double_precision(IIN,IGNORE_JUNK,xmax)
call read_value_integer(IIN,IGNORE_JUNK,nx)
call read_value_integer(IIN,IGNORE_JUNK,ngnod)
+ if ( ngnod == 9 .and. read_external_mesh ) then
+ print *, 'Number of control nodes must be equal to four when reading from external mesh.'
+ print *, 'ngnod = 9 is not yet supported.'
+ stop
+ end if
+
call read_value_logical(IIN,IGNORE_JUNK,initialfield)
call read_value_logical(IIN,IGNORE_JUNK,assign_external_model)
call read_value_logical(IIN,IGNORE_JUNK,TURN_ANISOTROPY_ON)
call read_value_logical(IIN,IGNORE_JUNK,TURN_ATTENUATION_ON)
+
+ print *, 'POYOP'
+ if ( read_external_mesh ) then
+ call read_mesh(mesh_file, nelmnts, elmnts, nnodes, num_start)
+
+ else
+ ! get interface data from external file to count the spectral elements along Z
+ print *,'Reading interface data from file DATA/',interfacesfile(1:len_trim(interfacesfile)),' to count the spectral elements'
+ open(unit=IIN_INTERFACES,file='DATA/'//interfacesfile,status='old')
+
+ max_npoints_interface = -1
+
+ ! read number of interfaces
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
+ if(number_of_interfaces < 2) stop 'not enough interfaces (minimum is 2)'
-! get interface data from external file to count the spectral elements along Z
- print *,'Reading interface data from file DATA/',interfacesfile(1:len_trim(interfacesfile)),' to count the spectral elements'
- open(unit=IIN_INTERFACES,file='DATA/'//interfacesfile,status='old')
+ ! loop on all the interfaces
+ do interface_current = 1,number_of_interfaces
+
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
+ if(npoints_interface_bottom < 2) stop 'not enough interface points (minimum is 2)'
+ max_npoints_interface = max(npoints_interface_bottom,max_npoints_interface)
+ print *,'Reading ',npoints_interface_bottom,' points for interface ',interface_current
+
+ ! loop on all the points describing this interface
+ do ipoint_current = 1,npoints_interface_bottom
+ call read_two_interface_points(IIN_INTERFACES,DONT_IGNORE_JUNK,xinterface_dummy,zinterface_dummy)
+ if(ipoint_current > 1 .and. xinterface_dummy <= xinterface_dummy_previous) &
+ stop 'interface points must be sorted in increasing X'
+ xinterface_dummy_previous = xinterface_dummy
+ enddo
+
+ enddo
+
+ ! define number of layers
+ number_of_layers = number_of_interfaces - 1
+
+ allocate(nz_layer(number_of_layers))
+
+ ! loop on all the layers
+ do ilayer = 1,number_of_layers
+
+ ! read number of spectral elements in vertical direction in this layer
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,nz_layer(ilayer))
+ if(nz_layer(ilayer) < 1) stop 'not enough spectral elements along Z in layer (minimum is 1)'
+ print *,'There are ',nz_layer(ilayer),' spectral elements along Z in layer ',ilayer
+
+ enddo
- max_npoints_interface = -1
+ close(IIN_INTERFACES)
+
+ ! compute total number of spectral elements in vertical direction
+ nz = sum(nz_layer)
+
+ print *
+ print *,'Total number of spectral elements along Z = ',nz
+ print *
+
+ nxread = nx
+ nzread = nz
+
+ ! multiply by 2 if elements have 9 nodes
+ if(ngnod == 9) then
+ nx = nx * 2
+ nz = nz * 2
+ nz_layer(:) = nz_layer(:) * 2
+ endif
+
+ nelmnts = nxread * nzread
+ allocate(elmnts(0:ngnod*nelmnts-1))
+ if ( ngnod == 4 ) then
+ num_elmnt = 0
+ do j = 1, nzread
+ do i = 1, nxread
+ elmnts(num_elmnt*ngnod) = (j-1)*(nxread+1) + (i-1)
+ elmnts(num_elmnt*ngnod+1) = (j-1)*(nxread+1) + (i-1) + 1
+ elmnts(num_elmnt*ngnod+2) = j*(nxread+1) + (i-1) + 1
+ elmnts(num_elmnt*ngnod+3) = j*(nxread+1) + (i-1)
+ num_elmnt = num_elmnt + 1
+ end do
+ end do
+ print *, 'POY', minval(elmnts)
+ else
+ num_elmnt = 0
+ do j = 1, nzread
+ do i = 1, nxread
+ elmnts(num_elmnt*ngnod) = (j-1)*(nxread+1) + (i-1)
+ elmnts(num_elmnt*ngnod+1) = (j-1)*(nxread+1) + (i-1) + 1
+ elmnts(num_elmnt*ngnod+2) = j*(nxread+1) + (i-1) + 1
+ elmnts(num_elmnt*ngnod+3) = j*(nxread+1) + (i-1)
+ elmnts(num_elmnt*ngnod+4) = (nxread+1)*(nzread+1) + (j-1)*nxread + (i-1)
+ elmnts(num_elmnt*ngnod+5) = (nxread+1)*(nzread+1) + nxread*(nzread+1) + (j-1)*(nxread*2+1) + (i-1)*2 + 2
+ elmnts(num_elmnt*ngnod+6) = (nxread+1)*(nzread+1) + j*nxread + (i-1)
+ elmnts(num_elmnt*ngnod+7) = (nxread+1)*(nzread+1) + nxread*(nzread+1) + (j-1)*(nxread*2+1) + (i-1)*2
+ elmnts(num_elmnt*ngnod+8) = (nxread+1)*(nzread+1) + nxread*(nzread+1) + (j-1)*(nxread*2+1) + (i-1)*2 + 1
+ num_elmnt = num_elmnt + 1
+ end do
+ end do
+
+ end if
+ !nnodes = (nx+1) * (nz+1)
+ end if
-! read number of interfaces
- call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
- if(number_of_interfaces < 2) stop 'not enough interfaces (minimum is 2)'
-
-! loop on all the interfaces
- do interface_current = 1,number_of_interfaces
-
- call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
- if(npoints_interface_bottom < 2) stop 'not enough interface points (minimum is 2)'
- max_npoints_interface = max(npoints_interface_bottom,max_npoints_interface)
- print *,'Reading ',npoints_interface_bottom,' points for interface ',interface_current
-
-! loop on all the points describing this interface
- do ipoint_current = 1,npoints_interface_bottom
- call read_two_interface_points(IIN_INTERFACES,DONT_IGNORE_JUNK,xinterface_dummy,zinterface_dummy)
- if(ipoint_current > 1 .and. xinterface_dummy <= xinterface_dummy_previous) &
- stop 'interface points must be sorted in increasing X'
- xinterface_dummy_previous = xinterface_dummy
- enddo
-
- enddo
-
-! define number of layers
- number_of_layers = number_of_interfaces - 1
-
- allocate(nz_layer(number_of_layers))
-
-! loop on all the layers
- do ilayer = 1,number_of_layers
-
-! read number of spectral elements in vertical direction in this layer
- call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,nz_layer(ilayer))
- if(nz_layer(ilayer) < 1) stop 'not enough spectral elements along Z in layer (minimum is 1)'
- print *,'There are ',nz_layer(ilayer),' spectral elements along Z in layer ',ilayer
-
- enddo
-
- close(IIN_INTERFACES)
-
-! compute total number of spectral elements in vertical direction
- nz = sum(nz_layer)
-
- print *
- print *,'Total number of spectral elements along Z = ',nz
- print *
-
- nxread = nx
- nzread = nz
-
-! multiply by 2 if elements have 9 nodes
- if(ngnod == 9) then
- nx = nx * 2
- nz = nz * 2
- nz_layer(:) = nz_layer(:) * 2
- endif
-
! read absorbing boundaries parameters
+ call read_value_logical(IIN,IGNORE_JUNK,any_abs)
call read_value_logical(IIN,IGNORE_JUNK,absbottom)
call read_value_logical(IIN,IGNORE_JUNK,absright)
call read_value_logical(IIN,IGNORE_JUNK,abstop)
call read_value_logical(IIN,IGNORE_JUNK,absleft)
+ if ( .not. any_abs ) then
+ absbottom = .false.
+ absright = .false.
+ abstop = .false.
+ absleft = .false.
+ end if
! read time step parameters
call read_value_integer(IIN,IGNORE_JUNK,nt)
@@ -302,7 +465,7 @@
allocate(cs(nb_materials))
allocate(aniso3(nb_materials))
allocate(aniso4(nb_materials))
- allocate(num_material(nx,nz))
+ allocate(num_material(nelmnts))
icodemat(:) = 0
rho(:) = 0.d0
@@ -310,7 +473,7 @@
cs(:) = 0.d0
aniso3(:) = 0.d0
aniso4(:) = 0.d0
- num_material(:,:) = 0
+ num_material(:) = 0
do imaterial=1,nb_materials
call read_material_parameters(IIN,DONT_IGNORE_JUNK,i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read)
@@ -325,7 +488,7 @@
aniso3(i) = aniso3read
aniso4(i) = aniso4read
enddo
-
+
print *
print *, 'Nb of solid or fluid materials = ',nb_materials
print *
@@ -345,71 +508,70 @@
print *
enddo
-! read the material numbers for each region
- call read_value_integer(IIN,IGNORE_JUNK,nbregion)
+
+ if ( read_external_mesh ) then
+ call read_mat(materials_file, nelmnts, num_material)
+ else
+ ! read the material numbers for each region
+ call read_value_integer(IIN,IGNORE_JUNK,nbregion)
- if(nbregion <= 0) stop 'Negative number of regions not allowed!'
+ if(nbregion <= 0) stop 'Negative number of regions not allowed!'
- print *
- print *, 'Nb of regions in the mesh = ',nbregion
- print *
+ print *
+ print *, 'Nb of regions in the mesh = ',nbregion
+ print *
- do iregion = 1,nbregion
+ do iregion = 1,nbregion
- call read_region_coordinates(IIN,DONT_IGNORE_JUNK,ixdebregion,ixfinregion,izdebregion,izfinregion,imaterial_number)
+ call read_region_coordinates(IIN,DONT_IGNORE_JUNK,ixdebregion,ixfinregion,izdebregion,izfinregion,imaterial_number)
- if(imaterial_number < 1) stop 'Negative material number not allowed!'
- if(ixdebregion < 1) stop 'Left coordinate of region negative!'
- if(ixfinregion > nxread) stop 'Right coordinate of region too high!'
- if(izdebregion < 1) stop 'Bottom coordinate of region negative!'
- if(izfinregion > nzread) stop 'Top coordinate of region too high!'
+ if(imaterial_number < 1) stop 'Negative material number not allowed!'
+ if(ixdebregion < 1) stop 'Left coordinate of region negative!'
+ if(ixfinregion > nxread) stop 'Right coordinate of region too high!'
+ if(izdebregion < 1) stop 'Bottom coordinate of region negative!'
+ if(izfinregion > nzread) stop 'Top coordinate of region too high!'
- print *,'Region ',iregion
- print *,'IX from ',ixdebregion,' to ',ixfinregion
- print *,'IZ from ',izdebregion,' to ',izfinregion
+ print *,'Region ',iregion
+ print *,'IX from ',ixdebregion,' to ',ixfinregion
+ print *,'IZ from ',izdebregion,' to ',izfinregion
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL) then
- vpregion = cp(imaterial_number)
- vsregion = cs(imaterial_number)
- print *,'Material # ',imaterial_number,' isotropic'
- if(vsregion < TINYVAL) then
- print *,'Material is fluid'
- else
- print *,'Material is solid'
- endif
- print *,'vp = ',vpregion
- print *,'vs = ',vsregion
- print *,'rho = ',rho(imaterial_number)
- poisson_ratio = 0.5d0*(vpregion*vpregion-2.d0*vsregion*vsregion) / (vpregion*vpregion-vsregion*vsregion)
- print *,'Poisson''s ratio = ',poisson_ratio
- if(poisson_ratio <= -1.00001d0 .or. poisson_ratio >= 0.50001d0) stop 'incorrect value of Poisson''s ratio'
- else
- print *,'Material # ',imaterial_number,' anisotropic'
- print *,'c11 = ',cp(imaterial_number)
- print *,'c13 = ',cs(imaterial_number)
- print *,'c33 = ',aniso3(imaterial_number)
- print *,'c44 = ',aniso4(imaterial_number)
- print *,'rho = ',rho(imaterial_number)
- endif
- print *,' -----'
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL) then
+ vpregion = cp(imaterial_number)
+ vsregion = cs(imaterial_number)
+ print *,'Material # ',imaterial_number,' isotropic'
+ if(vsregion < TINYVAL) then
+ print *,'Material is fluid'
+ else
+ print *,'Material is solid'
+ endif
+ print *,'vp = ',vpregion
+ print *,'vs = ',vsregion
+ print *,'rho = ',rho(imaterial_number)
+ poisson_ratio = 0.5d0*(vpregion*vpregion-2.d0*vsregion*vsregion) / (vpregion*vpregion-vsregion*vsregion)
+ print *,'Poisson''s ratio = ',poisson_ratio
+ if(poisson_ratio <= -1.00001d0 .or. poisson_ratio >= 0.50001d0) stop 'incorrect value of Poisson''s ratio'
+ else
+ print *,'Material # ',imaterial_number,' anisotropic'
+ print *,'c11 = ',cp(imaterial_number)
+ print *,'c13 = ',cs(imaterial_number)
+ print *,'c33 = ',aniso3(imaterial_number)
+ print *,'c44 = ',aniso4(imaterial_number)
+ print *,'rho = ',rho(imaterial_number)
+ endif
+ print *,' -----'
-! store density and velocity model
- do i = ixdebregion,ixfinregion
- do j = izdebregion,izfinregion
- if(ngnod == 4) then
- num_material(i,j) = imaterial_number
- else
- num_material(2*(i-1)+1,2*(j-1)+1) = imaterial_number
- num_material(2*(i-1)+1,2*(j-1)+2) = imaterial_number
- num_material(2*(i-1)+2,2*(j-1)+1) = imaterial_number
- num_material(2*(i-1)+2,2*(j-1)+2) = imaterial_number
- endif
+ ! store density and velocity model
+ do i = ixdebregion,ixfinregion
+ do j = izdebregion,izfinregion
+ num_material((j-1)*nxread+i) = imaterial_number
+ enddo
+ enddo
+
enddo
- enddo
- enddo
+ if(minval(num_material) <= 0) stop 'Velocity model not entirely set...'
- if(minval(num_material) <= 0) stop 'Velocity model not entirely set...'
+ end if
close(IIN)
@@ -421,369 +583,659 @@
if(ngnod /= 4 .and. ngnod /= 9) stop 'ngnod different from 4 or 9!'
print *
- if(ngnod == 4) then
- print *,'The mesh contains ',nx,' x ',nz,' elements'
- else
- print *,'The mesh contains ',nx/2,' x ',nz/2,' elements'
- endif
+ print *,'The mesh contains ',nelmnts,' elements'
print *
print *,'Control elements have ',ngnod,' nodes'
print *
!---
-! allocate arrays for the grid
- allocate(x(0:nx,0:nz))
- allocate(z(0:nx,0:nz))
+ if ( .not. read_external_mesh ) then
+ ! allocate arrays for the grid
+ allocate(x(0:nx,0:nz))
+ allocate(z(0:nx,0:nz))
- x(:,:) = 0.d0
- z(:,:) = 0.d0
+ x(:,:) = 0.d0
+ z(:,:) = 0.d0
-! get interface data from external file
- print *,'Reading interface data from file DATA/',interfacesfile(1:len_trim(interfacesfile))
- open(unit=IIN_INTERFACES,file='DATA/'//interfacesfile,status='old')
+ ! get interface data from external file
+ print *,'Reading interface data from file DATA/',interfacesfile(1:len_trim(interfacesfile))
+ open(unit=IIN_INTERFACES,file='DATA/'//interfacesfile,status='old')
- allocate(xinterface_bottom(max_npoints_interface))
- allocate(zinterface_bottom(max_npoints_interface))
- allocate(coefs_interface_bottom(max_npoints_interface))
+ allocate(xinterface_bottom(max_npoints_interface))
+ allocate(zinterface_bottom(max_npoints_interface))
+ allocate(coefs_interface_bottom(max_npoints_interface))
- allocate(xinterface_top(max_npoints_interface))
- allocate(zinterface_top(max_npoints_interface))
- allocate(coefs_interface_top(max_npoints_interface))
+ allocate(xinterface_top(max_npoints_interface))
+ allocate(zinterface_top(max_npoints_interface))
+ allocate(coefs_interface_top(max_npoints_interface))
-! read number of interfaces
- call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
+ ! read number of interfaces
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
-! read bottom interface
- call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
+ ! read bottom interface
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
-! loop on all the points describing this interface
- do ipoint_current = 1,npoints_interface_bottom
- call read_two_interface_points(IIN_INTERFACES,DONT_IGNORE_JUNK, &
+ ! loop on all the points describing this interface
+ do ipoint_current = 1,npoints_interface_bottom
+ call read_two_interface_points(IIN_INTERFACES,DONT_IGNORE_JUNK, &
xinterface_bottom(ipoint_current),zinterface_bottom(ipoint_current))
- enddo
+ enddo
-! loop on all the layers
- do ilayer = 1,number_of_layers
+ ! loop on all the layers
+ do ilayer = 1,number_of_layers
-! read top interface
- call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_top)
+ ! read top interface
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_top)
-! loop on all the points describing this interface
- do ipoint_current = 1,npoints_interface_top
- call read_two_interface_points(IIN_INTERFACES,DONT_IGNORE_JUNK, &
- xinterface_top(ipoint_current),zinterface_top(ipoint_current))
- enddo
+ ! loop on all the points describing this interface
+ do ipoint_current = 1,npoints_interface_top
+ call read_two_interface_points(IIN_INTERFACES,DONT_IGNORE_JUNK, &
+ xinterface_top(ipoint_current),zinterface_top(ipoint_current))
+ enddo
-! compute the spline for the bottom interface, impose the tangent on both edges
- tang1 = (zinterface_bottom(2)-zinterface_bottom(1)) / (xinterface_bottom(2)-xinterface_bottom(1))
- tangN = (zinterface_bottom(npoints_interface_bottom)-zinterface_bottom(npoints_interface_bottom-1)) / &
- (xinterface_bottom(npoints_interface_bottom)-xinterface_bottom(npoints_interface_bottom-1))
- call spline(xinterface_bottom,zinterface_bottom,npoints_interface_bottom,tang1,tangN,coefs_interface_bottom)
+ ! compute the spline for the bottom interface, impose the tangent on both edges
+ tang1 = (zinterface_bottom(2)-zinterface_bottom(1)) / (xinterface_bottom(2)-xinterface_bottom(1))
+ tangN = (zinterface_bottom(npoints_interface_bottom)-zinterface_bottom(npoints_interface_bottom-1)) / &
+ (xinterface_bottom(npoints_interface_bottom)-xinterface_bottom(npoints_interface_bottom-1))
+ call spline(xinterface_bottom,zinterface_bottom,npoints_interface_bottom,tang1,tangN,coefs_interface_bottom)
-! compute the spline for the top interface, impose the tangent on both edges
- tang1 = (zinterface_top(2)-zinterface_top(1)) / (xinterface_top(2)-xinterface_top(1))
- tangN = (zinterface_top(npoints_interface_top)-zinterface_top(npoints_interface_top-1)) / &
- (xinterface_top(npoints_interface_top)-xinterface_top(npoints_interface_top-1))
- call spline(xinterface_top,zinterface_top,npoints_interface_top,tang1,tangN,coefs_interface_top)
+ ! compute the spline for the top interface, impose the tangent on both edges
+ tang1 = (zinterface_top(2)-zinterface_top(1)) / (xinterface_top(2)-xinterface_top(1))
+ tangN = (zinterface_top(npoints_interface_top)-zinterface_top(npoints_interface_top-1)) / &
+ (xinterface_top(npoints_interface_top)-xinterface_top(npoints_interface_top-1))
+ call spline(xinterface_top,zinterface_top,npoints_interface_top,tang1,tangN,coefs_interface_top)
-! check if we are in the last layer, which contains topography,
-! and modify the position of the source accordingly if it is located exactly at the surface
- if(source_surf .and. ilayer == number_of_layers) &
- zs = value_spline(xs,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
+ ! check if we are in the last layer, which contains topography,
+ ! and modify the position of the source accordingly if it is located exactly at the surface
+ if(source_surf .and. ilayer == number_of_layers) &
+ zs = value_spline(xs,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
-! compute the offset of this layer in terms of number of spectral elements below along Z
- if(ilayer > 1) then
- ioffset = sum(nz_layer(1:ilayer-1))
- else
- ioffset = 0
- endif
+ ! compute the offset of this layer in terms of number of spectral elements below along Z
+ if(ilayer > 1) then
+ ioffset = sum(nz_layer(1:ilayer-1))
+ else
+ ioffset = 0
+ endif
-!--- definition of the mesh
+ !--- definition of the mesh
- do ix = 0,nx
+ do ix = 0,nx
-! evenly spaced points along X
- absx = xmin + (xmax - xmin) * dble(ix) / dble(nx)
+ ! evenly spaced points along X
+ absx = xmin + (xmax - xmin) * dble(ix) / dble(nx)
-! value of the bottom and top splines
- bot0 = value_spline(absx,xinterface_bottom,zinterface_bottom,coefs_interface_bottom,npoints_interface_bottom)
- top0 = value_spline(absx,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
+ ! value of the bottom and top splines
+ bot0 = value_spline(absx,xinterface_bottom,zinterface_bottom,coefs_interface_bottom,npoints_interface_bottom)
+ top0 = value_spline(absx,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
- do iz = 0,nz_layer(ilayer)
+ do iz = 0,nz_layer(ilayer)
-! linear interpolation between bottom and top
- gamma = dble(iz) / dble(nz_layer(ilayer))
- a00 = 1.d0 - gamma
- a01 = gamma
+ ! linear interpolation between bottom and top
+ gamma = dble(iz) / dble(nz_layer(ilayer))
+ a00 = 1.d0 - gamma
+ a01 = gamma
-! coordinates of the grid points
- x(ix,iz + ioffset) = absx
- z(ix,iz + ioffset) = a00*bot0 + a01*top0
+ ! coordinates of the grid points
+ x(ix,iz + ioffset) = absx
+ z(ix,iz + ioffset) = a00*bot0 + a01*top0
- enddo
+ enddo
- enddo
+ enddo
-! the top interface becomes the bottom interface before switching to the next layer
- npoints_interface_bottom = npoints_interface_top
- xinterface_bottom(:) = xinterface_top(:)
- zinterface_bottom(:) = zinterface_top(:)
+ ! the top interface becomes the bottom interface before switching to the next layer
+ npoints_interface_bottom = npoints_interface_top
+ xinterface_bottom(:) = xinterface_top(:)
+ zinterface_bottom(:) = zinterface_top(:)
- enddo
+ enddo
- close(IIN_INTERFACES)
+ close(IIN_INTERFACES)
+
+ nnodes = (nz+1)*(nx+1)
+ allocate(nodes_coords(2,nnodes))
+ if ( ngnod == 4 ) then
+ do j = 0, nz
+ do i = 0, nx
+ num_node = num_4(i,j,nxread)
+ nodes_coords(1, num_node) = x(i,j)
+ nodes_coords(2, num_node) = z(i,j)
+
+ end do
+ end do
+
+ else
+ do j = 0, nz
+ do i = 0, nx
+ num_node = num_9(i,j,nxread,nzread)
+ nodes_coords(1, num_node) = x(i,j)
+ nodes_coords(2, num_node) = z(i,j)
+
+ end do
+ end do
+
+ end if
+ else
+ call read_nodes_coords(nodes_coords_file, nnodes, nodes_coords)
+ end if
+
+ if ( read_external_mesh ) then
+ call read_acoustic_surface(free_surface_file, nelem_acoustic_surface, acoustic_surface, &
+ nelmnts, num_material, ANISOTROPIC_MATERIAL, nb_materials, icodemat, cs, num_start)
+
+ if ( any_abs ) then
+ call read_abs_surface(absorbing_surface_file, nelemabs, abs_surface, num_start)
+ end if
+
+ else
+ ! count the number of acoustic free-surface elements
+ nelem_acoustic_surface = 0
+ j = nzread
+ do i = 1,nxread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ end if
+ enddo
+
+ allocate(acoustic_surface(4,nelem_acoustic_surface))
+
+ nelem_acoustic_surface = 0
+ j = nzread
+ do i = 1,nxread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+ acoustic_surface(2,nelem_acoustic_surface) = 2
+ acoustic_surface(3,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
+ acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
+ end if
+ end do
+
+ !
+ !--- definition of absorbing boundaries
+ !
+ nelemabs = 0
+ if(absbottom) nelemabs = nelemabs + nxread
+ if(abstop) nelemabs = nelemabs + nxread
+ if(absleft) nelemabs = nelemabs + nzread
+ if(absright) nelemabs = nelemabs + nzread
+
+ allocate(abs_surface(4,nelemabs))
+
+ ! generate the list of absorbing elements
+ if(nelemabs > 0) then
+ nelemabs = 0
+ do iz = 1,nzread
+ do ix = 1,nxread
+ inumelem = (iz-1)*nxread + ix
+ if(absbottom .and. iz == 1) then
+ nelemabs = nelemabs + 1
+ abs_surface(1,nelemabs) = inumelem-1
+ abs_surface(2,nelemabs) = 2
+ abs_surface(3,nelemabs) = elmnts(0+ngnod*(inumelem-1))
+ abs_surface(4,nelemabs) = elmnts(1+ngnod*(inumelem-1))
+ end if
+ if(absright .and. ix == nxread) then
+ nelemabs = nelemabs + 1
+ abs_surface(1,nelemabs) = inumelem-1
+ abs_surface(2,nelemabs) = 2
+ abs_surface(3,nelemabs) = elmnts(1+ngnod*(inumelem-1))
+ abs_surface(4,nelemabs) = elmnts(2+ngnod*(inumelem-1))
+ end if
+ if(abstop .and. iz == nzread) then
+ nelemabs = nelemabs + 1
+ abs_surface(1,nelemabs) = inumelem-1
+ abs_surface(2,nelemabs) = 2
+ abs_surface(3,nelemabs) = elmnts(3+ngnod*(inumelem-1))
+ abs_surface(4,nelemabs) = elmnts(2+ngnod*(inumelem-1))
+ end if
+ if(absleft .and. ix == 1) then
+ nelemabs = nelemabs + 1
+ abs_surface(1,nelemabs) = inumelem-1
+ abs_surface(2,nelemabs) = 2
+ abs_surface(3,nelemabs) = elmnts(0+ngnod*(inumelem-1))
+ abs_surface(4,nelemabs) = elmnts(3+ngnod*(inumelem-1))
+ end if
+ end do
+ end do
+ end if
+
+ end if
+
+
! compute min and max of X and Z in the grid
print *
- print *,'Min and max value of X in the grid = ',minval(x),maxval(x)
- print *,'Min and max value of Z in the grid = ',minval(z),maxval(z)
+ print *,'Min and max value of X in the grid = ',minval(nodes_coords(1,:)),maxval(nodes_coords(1,:))
+ print *,'Min and max value of Z in the grid = ',minval(nodes_coords(2,:)),maxval(nodes_coords(2,:))
print *
-
+
+
! ***
! *** create a Gnuplot file that displays the grid
! ***
- print *
- print *,'Saving the grid in Gnuplot format...'
+!!$ print *
+!!$ print *,'Saving the grid in Gnuplot format...'
+!!$
+!!$ open(unit=20,file='OUTPUT_FILES/gridfile.gnu',status='unknown')
+!!$
+!!$! draw horizontal lines of the grid
+!!$ print *,'drawing horizontal lines of the grid'
+!!$ istepx = 1
+!!$ if(ngnod == 4) then
+!!$ istepz = 1
+!!$ else
+!!$ istepz = 2
+!!$ endif
+!!$ do ili=0,nz,istepz
+!!$ do icol=0,nx-istepx,istepx
+!!$ write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
+!!$ write(20,*) sngl(x(icol+istepx,ili)),sngl(z(icol+istepx,ili))
+!!$ write(20,10)
+!!$ enddo
+!!$ enddo
+!!$
+!!$! draw vertical lines of the grid
+!!$ print *,'drawing vertical lines of the grid'
+!!$ if(ngnod == 4) then
+!!$ istepx = 1
+!!$ else
+!!$ istepx = 2
+!!$ endif
+!!$ istepz = 1
+!!$ do icol=0,nx,istepx
+!!$ do ili=0,nz-istepz,istepz
+!!$ write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
+!!$ write(20,*) sngl(x(icol,ili+istepz)),sngl(z(icol,ili+istepz))
+!!$ write(20,10)
+!!$ enddo
+!!$ enddo
+!!$
+!!$ 10 format('')
+!!$
+!!$ close(20)
+!!$
+!!$! create a Gnuplot script to display the grid
+!!$ open(unit=20,file='OUTPUT_FILES/plotgnu',status='unknown')
+!!$ write(20,*) '#set term X11'
+!!$ write(20,*) 'set term postscript landscape monochrome solid "Helvetica" 22'
+!!$ write(20,*) 'set output "grid.ps"'
+!!$ write(20,*) '#set xrange [',sngl(minval(x)),':',sngl(maxval(x)),']'
+!!$ write(20,*) '#set yrange [',sngl(minval(z)),':',sngl(maxval(z)),']'
+!!$! use same unit length on both X and Y axes
+!!$ write(20,*) 'set size ratio -1'
+!!$ write(20,*) 'plot "gridfile.gnu" title "Macrobloc mesh" w l'
+!!$ write(20,*) 'pause -1 "Hit any key..."'
+!!$ close(20)
+!!$
+!!$ print *,'Grid saved in Gnuplot format...'
+!!$ print *
+
+
+
+ !*****************************
+ ! Partitionning
+ !*****************************
+ allocate(part(0:nelmnts-1))
- open(unit=20,file='OUTPUT_FILES/gridfile.gnu',status='unknown')
-
-! draw horizontal lines of the grid
- print *,'drawing horizontal lines of the grid'
- istepx = 1
- if(ngnod == 4) then
- istepz = 1
+!!$ if ( nproc == 1 ) then
+!!$ ! There is only one process; no need for partitionning
+!!$ call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
+!!$ part(:) = 0
+!!$
+!!$ else
+ if ( ngnod == 9 ) then
+ allocate(elmnts_bis(0:ESIZE*nelmnts-1))
+ do i = 0, nelmnts-1
+ elmnts_bis(i*esize:i*esize+esize-1) = elmnts(i*ngnod:i*ngnod+esize-1)
+ end do
+
+ call mesh2dual_ncommonnodes(nelmnts, (nxread+1)*(nzread+1), elmnts_bis, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
+
else
- istepz = 2
- endif
- do ili=0,nz,istepz
- do icol=0,nx-istepx,istepx
- write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
- write(20,*) sngl(x(icol+istepx,ili)),sngl(z(icol+istepx,ili))
- write(20,10)
- enddo
- enddo
+ call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
+
+ end if
+
+ nb_edges = xadj(nelmnts)
+
+ call read_weights(nelmnts, vwgt, nb_edges, adjwgt)
+
+ if ( nproc == 1 ) then
+ part(:) = 0
+ else
-! draw vertical lines of the grid
- print *,'drawing vertical lines of the grid'
- if(ngnod == 4) then
- istepx = 1
+ select case (partitionning_method)
+ case(1)
+ do iproc = 0, nproc-2
+ part(iproc*floor(real(nelmnts)/real(nproc)):(iproc+1)*floor(real(nelmnts)/real(nproc))-1) = iproc
+ end do
+ part(floor(real(nelmnts)/real(nproc))*(nproc-1):nelmnts-1) = nproc - 1
+!!$ part(0:1659) = 0
+!!$ part(1660:3258) = 1
+!!$ part(3259:4799) = 2
+!!$ print *, 'WWWWWW', nelmnts-1
+
+ case(2)
+#ifdef USE_METIS
+ call Part_metis(nelmnts, xadj, adjncy, vwgt, adjwgt, nproc, nb_edges, edgecut, part, metis_options)
+#else
+ print *, 'This version of SPECFEM was not compiled with support of METIS.'
+ print *, 'Please recompile with -DUSE_METIS in order to enable use of METIS.'
+ stop
+#endif
+
+ case(3)
+#ifdef USE_SCOTCH
+ call Part_scotch(nelmnts, xadj, adjncy, vwgt, adjwgt, nproc, nb_edges, edgecut, part, scotch_strategy)
+#else
+ print *, 'This version of SPECFEM was not compiled with support of SCOTCH.'
+ print *, 'Please recompile with -DUSE_SCOTCH in order to enable use of SCOTCH.'
+ stop
+#endif
+
+ end select
+
+ end if
+
+ ! beware of fluid solid edges
+ if ( ngnod == 9 ) then
+ call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts_bis, nb_materials, cs, num_material, &
+ nproc, part, nedges_coupled, edges_coupled)
else
- istepx = 2
- endif
- istepz = 1
- do icol=0,nx,istepx
- do ili=0,nz-istepz,istepz
- write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
- write(20,*) sngl(x(icol,ili+istepz)),sngl(z(icol,ili+istepz))
- write(20,10)
- enddo
- enddo
+ call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, nb_materials, cs, num_material, &
+ nproc, part, nedges_coupled, edges_coupled)
+ end if
+
+ call Construct_glob2loc_elmnts(nelmnts, part, nproc, glob2loc_elmnts)
+
+ if ( ngnod == 9 ) then
+!!$ print *, 'POUMPOUM', nnodes, elmnts(0:8)
+!!$ print *, nodes_coords(1,elmnts(0:8)+1)
+!!$ print *, nodes_coords(2,elmnts(0:8)+1)
+!!$ print *, x(1,0), z(1,0)
+!!$ print *, x(1,1), z(1,1)
+!!$ print *, num_9(2,1,nxread,nzread), x(2,1), z(2,1)
+!!$ print *, num_9(0,1,nxread,nzread),x(0,1), z(0,1)
+
+ deallocate(nnodes_elmnts)
+ deallocate(nodes_elmnts)
+ allocate(nnodes_elmnts(0:nnodes-1))
+ allocate(nodes_elmnts(0:nsize*nnodes-1))
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+ do i = 0, ngnod*nelmnts-1
+ nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
+ nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
+
+ end do
+ end if
+
+ call Construct_glob2loc_nodes(nelmnts, nnodes, nnodes_elmnts, nodes_elmnts, part, nproc, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
- 10 format('')
+ if ( nproc /= 1 ) then
+ if ( ngnod == 9 ) then
+ call Construct_interfaces(nelmnts, nproc, part, elmnts_bis, xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, nb_materials, cs, num_material)
+ else
+ call Construct_interfaces(nelmnts, nproc, part, elmnts, xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, nb_materials, cs, num_material)
+ end if
+ print *, '04'
+ allocate(my_interfaces(0:ninterfaces-1))
+ allocate(my_nb_interfaces(0:ninterfaces-1))
+ print *, '05'
+ end if
+
+ if ( any_abs ) then
+ call merge_abs_boundaries(nelemabs, nelemabs_merge, abs_surface, abs_surface_char, abs_surface_merge, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ nodes_coords, &
+ nedges_coupled, edges_coupled, nb_materials, cs, num_material, &
+ nelmnts, &
+ elmnts, ngnod)
+ print *, 'nelemabs_merge', nelemabs_merge
+ end if
- close(20)
-
-! create a Gnuplot script to display the grid
- open(unit=20,file='OUTPUT_FILES/plotgnu',status='unknown')
- write(20,*) '#set term X11'
- write(20,*) 'set term postscript landscape monochrome solid "Helvetica" 22'
- write(20,*) 'set output "grid.ps"'
- write(20,*) '#set xrange [',sngl(minval(x)),':',sngl(maxval(x)),']'
- write(20,*) '#set yrange [',sngl(minval(z)),':',sngl(maxval(z)),']'
-! use same unit length on both X and Y axes
- write(20,*) 'set size ratio -1'
- write(20,*) 'plot "gridfile.gnu" title "Macrobloc mesh" w l'
- write(20,*) 'pause -1 "Hit any key..."'
- close(20)
-
- print *,'Grid saved in Gnuplot format...'
- print *
-
! *** generate the database for the solver
- open(unit=15,file='OUTPUT_FILES/Database',status='unknown')
+ do iproc = 0, nproc-1
+
+ write(prname, "('/Database',i5.5)") iproc
+ open(unit=15,file='./OUTPUT_FILES'//prname,status='unknown')
- write(15,*) '#'
- write(15,*) '# Database for SPECFEM2D'
- write(15,*) '# Dimitri Komatitsch, (c) University of Pau, France'
- write(15,*) '#'
+ write(15,*) '#'
+ write(15,*) '# Database for SPECFEM2D'
+ write(15,*) '# Dimitri Komatitsch, (c) University of Pau, France'
+ write(15,*) '#'
+
+ write(15,*) 'Title of the simulation'
+ write(15,"(a50)") title
+
- write(15,*) 'Title of the simulation'
- write(15,"(a50)") title
+ call write_glob2loc_nodes_database(15, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 1)
+
- npgeo = (nx+1)*(nz+1)
- if(ngnod == 4) then
- nspec = nx*nz
- else
- nspec = nx*nz/4
- endif
- write(15,*) 'npgeo'
- write(15,*) npgeo
+ call write_partition_database(15, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 1)
+
- write(15,*) 'gnuplot interpol'
- write(15,*) gnuplot,interpol
+ write(15,*) 'npgeo'
+ write(15,*) npgeo
+
+ write(15,*) 'gnuplot interpol'
+ write(15,*) gnuplot,interpol
+
+ write(15,*) 'NTSTEP_BETWEEN_OUTPUT_INFO'
+ write(15,*) NTSTEP_BETWEEN_OUTPUT_INFO
+
+ write(15,*) 'output_postscript_snapshot output_color_image colors numbers'
+ write(15,*) output_postscript_snapshot,output_color_image,' 1 0'
- write(15,*) 'NTSTEP_BETWEEN_OUTPUT_INFO'
- write(15,*) NTSTEP_BETWEEN_OUTPUT_INFO
+ write(15,*) 'meshvect modelvect boundvect cutsnaps subsamp sizemax_arrows'
+ write(15,*) meshvect,modelvect,boundvect,cutsnaps,subsamp,sizemax_arrows
+
+ write(15,*) 'anglerec'
+ write(15,*) anglerec
+
+ write(15,*) 'initialfield'
+ write(15,*) initialfield
- write(15,*) 'output_postscript_snapshot output_color_image colors numbers'
- write(15,*) output_postscript_snapshot,output_color_image,' 1 0'
+ write(15,*) 'seismotype imagetype'
+ write(15,*) seismotype,imagetype
- write(15,*) 'meshvect modelvect boundvect cutsnaps subsamp sizemax_arrows'
- write(15,*) meshvect,modelvect,boundvect,cutsnaps,subsamp,sizemax_arrows
+ write(15,*) 'assign_external_model outputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
+ write(15,*) assign_external_model,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+
+ write(15,*) 'nt deltat'
+ write(15,*) nt,deltat
- write(15,*) 'anglerec'
- write(15,*) anglerec
+ write(15,*) 'source'
+ write(15,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
+
+ write(15,*) 'Coordinates of macrobloc mesh (coorg):'
+
- write(15,*) 'initialfield'
- write(15,*) initialfield
+ call write_glob2loc_nodes_database(15, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 2)
+
- write(15,*) 'seismotype imagetype'
- write(15,*) seismotype,imagetype
+ write(15,*) 'numat ngnod nspec pointsdisp plot_lowerleft_corner_only'
+ write(15,*) nb_materials,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
- write(15,*) 'assign_external_model outputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
- write(15,*) assign_external_model,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+
+ !call Write_surface_database(15, nelemabs, abs_surface, nelemabs_loc, &
+ ! nproc, iproc, glob2loc_elmnts, &
+ ! glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 1)
+ if ( any_abs ) then
+ call write_abs_merge_database(15, nelemabs_merge, nelemabs_loc, &
+ abs_surface_char, abs_surface_merge, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ glob2loc_elmnts, part, iproc, 1)
+ else
+ nelemabs_loc = 0
+ end if
+
+ call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
+ nproc, iproc, glob2loc_elmnts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 1)
+
- write(15,*) 'nt deltat'
- write(15,*) nt,deltat
+ call write_acoustic_elastic_edges_database(15, nedges_coupled, nedges_coupled_loc, &
+ edges_coupled, glob2loc_elmnts, nelmnts, part, iproc, 1)
+
+ write(15,*) 'nelemabs nelem_acoustic_surface num_fluid_solid_edges'
+ write(15,*) nelemabs_loc,nelem_acoustic_surface_loc,nedges_coupled_loc
+
+
+ write(15,*) 'Material sets (num 1 rho vp vs 0 0) or (num 2 rho c11 c13 c33 c44)'
+ do i=1,nb_materials
+ write(15,*) i,icodemat(i),rho(i),cp(i),cs(i),aniso3(i),aniso4(i)
+ enddo
+
+ write(15,*) 'Arrays kmato and knods for each bloc:'
- write(15,*) 'source'
- write(15,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
+
+ call write_partition_database(15, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 2)
+
+ if ( nproc /= 1 ) then
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 1)
+
+ write(15,*) 'Interfaces:'
+ write(15,*) my_ninterface, maxval(my_nb_interfaces)
+
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 2)
- write(15,*) 'Coordinates of macrobloc mesh (coorg):'
- do j=0,nz
- do i=0,nx
- write(15,*) num(i,j,nx),x(i,j),z(i,j)
- enddo
- enddo
+ else
+ write(15,*) 'Interfaces:'
+ write(15,*) 0, 0
+
+ end if
+
-!
-!--- definition of absorbing boundaries
-!
- nelemabs = 0
- if(absbottom) nelemabs = nelemabs + nx
- if(abstop) nelemabs = nelemabs + nx
- if(absleft) nelemabs = nelemabs + nz
- if(absright) nelemabs = nelemabs + nz
+ write(15,*) 'List of absorbing elements (bottom right top left):'
+ !call Write_surface_database(15, nelemabs, abs_surface, nelemabs_loc, &
+ ! nproc, iproc, glob2loc_elmnts, &
+ ! glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 2)
+ if ( any_abs ) then
+ call write_abs_merge_database(15, nelemabs_merge, nelemabs_loc, &
+ abs_surface_char, abs_surface_merge, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ glob2loc_elmnts, part, iproc, 2)
+ end if
+
+ write(15,*) 'List of acoustic free-surface elements:'
+ call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
+ nproc, iproc, glob2loc_elmnts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 2)
+
-! we have counted the elements twice if they have nine nodes
- if(ngnod == 9) nelemabs = nelemabs / 2
+ write(15,*) 'List of acoustic elastic coupled edges:'
+ call write_acoustic_elastic_edges_database(15, nedges_coupled, nedges_coupled_loc, &
+ edges_coupled, glob2loc_elmnts, nelmnts, part, iproc, 2)
+ end do
+
+
+!!$
+!!$ open(unit=15,file='OUTPUT_FILES/Database',status='unknown')
+!!$
-! also remove the corner elements, which have been counted twice
- if(absbottom .and. absleft) nelemabs = nelemabs - 1
- if(absbottom .and. absright) nelemabs = nelemabs - 1
- if(abstop .and. absleft) nelemabs = nelemabs - 1
- if(abstop .and. absright) nelemabs = nelemabs - 1
+!!$
+!!$ npgeo = nnodes
+!!$ nspec = nelmnts
+!!$
-! count the number of acoustic free-surface elements
- nelem_acoustic_surface = 0
- do j = 1,nzread
- do i = 1,nxread
- if(ngnod == 4) then
- imaterial_number = num_material(i,j)
- else
- imaterial_number = num_material(2*(i-1)+1,2*(j-1)+1)
- endif
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL &
- .and. j == nzread) nelem_acoustic_surface = nelem_acoustic_surface + 1
- enddo
- enddo
+!!$ do j=0,nz
+!!$ do i=0,nx
+!!$ write(15,*) num(i,j,nx),x(i,j),z(i,j)
+!!$ enddo
+!!$ enddo
+!!$
+!!$
+!!$
+!!$
+!!$
- write(15,*) 'numat ngnod nspec pointsdisp plot_lowerleft_corner_only'
- write(15,*) nb_materials,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
- write(15,*) 'nelemabs nelem_acoustic_surface'
- write(15,*) nelemabs,nelem_acoustic_surface
+!!$
- write(15,*) 'Material sets (num 1 rho vp vs 0 0) or (num 2 rho c11 c13 c33 c44)'
- do i=1,nb_materials
- write(15,*) i,icodemat(i),rho(i),cp(i),cs(i),aniso3(i),aniso4(i)
- enddo
-
-
- write(15,*) 'Arrays kmato and knods for each bloc:'
-
- k = 0
- if(ngnod == 4) then
- do j=0,nz-1
- do i=0,nx-1
- k = k + 1
- imaterial_number = num_material(i+1,j+1)
- write(15,*) k,imaterial_number,num(i,j,nx),num(i+1,j,nx),num(i+1,j+1,nx),num(i,j+1,nx)
- enddo
- enddo
- else
- do j=0,nz-2,2
- do i=0,nx-2,2
- k = k + 1
- imaterial_number = num_material(i+1,j+1)
- write(15,*) k,imaterial_number,num(i,j,nx),num(i+2,j,nx),num(i+2,j+2,nx),num(i,j+2,nx), &
- num(i+1,j,nx),num(i+2,j+1,nx),num(i+1,j+2,nx),num(i,j+1,nx),num(i+1,j+1,nx)
- enddo
- enddo
- endif
-
-!
-!--- save absorbing boundaries
-!
- print *
- print *,'There is a total of ',nelemabs,' absorbing elements'
- print *
- print *,'Active absorbing boundaries:'
- print *
- print *,'Bottom = ',absbottom
- print *,'Right = ',absright
- print *,'Top = ',abstop
- print *,'Left = ',absleft
- print *
-
-! generate the list of absorbing elements
- if(nelemabs > 0) then
- write(15,*) 'List of absorbing elements (bottom right top left):'
- do iz = 1,nzread
- do ix = 1,nxread
- codebottom = .false.
- coderight = .false.
- codetop = .false.
- codeleft = .false.
- inumelem = (iz-1)*nxread + ix
- if(absbottom .and. iz == 1) codebottom = .true.
- if(absright .and. ix == nxread) coderight = .true.
- if(abstop .and. iz == nzread) codetop = .true.
- if(absleft .and. ix == 1) codeleft = .true.
- if(codebottom .or. coderight .or. codetop .or. codeleft) write(15,*) inumelem,codebottom,coderight,codetop,codeleft
- enddo
- enddo
- endif
-
-!
-!--- save acoustic free-surface elements
-!
- print *
- print *,'There is a total of ',nelem_acoustic_surface,' acoustic free-surface elements'
- print *
-
-! generate the list of acoustic free-surface elements
- if(nelem_acoustic_surface > 0) then
- write(15,*) 'List of acoustic free-surface elements:'
- do j = 1,nzread
- do i = 1,nxread
- inumelem = (j-1)*nxread + i
- if(ngnod == 4) then
- imaterial_number = num_material(i,j)
- else
- imaterial_number = num_material(2*(i-1)+1,2*(j-1)+1)
- endif
-! in this simple mesher, it is always the top edge that is at the free surface
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL .and. j == nzread) &
- write(15,*) inumelem,ITOP
- enddo
- enddo
- endif
-
- close(15)
-
+!!$
+!!$ k = 0
+!!$ if(ngnod == 4) then
+!!$ do j=0,nz-1
+!!$ do i=0,nx-1
+!!$ k = k + 1
+!!$ imaterial_number = num_material(i+1,j+1)
+!!$ write(15,*) k,imaterial_number,num(i,j,nx),num(i+1,j,nx),num(i+1,j+1,nx),num(i,j+1,nx)
+!!$ enddo
+!!$ enddo
+!!$ else
+!!$ do j=0,nz-2,2
+!!$ do i=0,nx-2,2
+!!$ k = k + 1
+!!$ imaterial_number = num_material(i+1,j+1)
+!!$ write(15,*) k,imaterial_number,num(i,j,nx),num(i+2,j,nx),num(i+2,j+2,nx),num(i,j+2,nx), &
+!!$ num(i+1,j,nx),num(i+2,j+1,nx),num(i+1,j+2,nx),num(i,j+1,nx),num(i+1,j+1,nx)
+!!$ enddo
+!!$ enddo
+!!$ endif
+!!$
+!!$!
+!!$!--- save absorbing boundaries
+!!$!
+!!$ print *
+!!$ print *,'There is a total of ',nelemabs,' absorbing elements'
+!!$ print *
+!!$ print *,'Active absorbing boundaries:'
+!!$ print *
+!!$ print *,'Bottom = ',absbottom
+!!$ print *,'Right = ',absright
+!!$ print *,'Top = ',abstop
+!!$ print *,'Left = ',absleft
+!!$ print *
+!!$
+!!$
+!!$
+!!$!
+!!$!--- save acoustic free-surface elements
+!!$!
+!!$ print *
+!!$ print *,'There is a total of ',nelem_acoustic_surface,' acoustic free-surface elements'
+!!$ print *
+!!$
+!!$! generate the list of acoustic free-surface elements
+!!$ if(nelem_acoustic_surface > 0) then
+!!$ write(15,*) 'List of acoustic free-surface elements:'
+!!$ do j = 1,nzread
+!!$ do i = 1,nxread
+!!$ inumelem = (j-1)*nxread + i
+!!$ if(ngnod == 4) then
+!!$ imaterial_number = num_material(i,j)
+!!$ else
+!!$ imaterial_number = num_material(2*(i-1)+1,2*(j-1)+1)
+!!$ endif
+!!$! in this simple mesher, it is always the top edge that is at the free surface
+!!$ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL .and. j == nzread) &
+!!$ write(15,*) inumelem,ITOP
+!!$ enddo
+!!$ enddo
+!!$ endif
+!!$
+!!$ close(15)
+!!$
! print position of the source
print *
print *,'Position (x,z) of the source = ',xs,zs
@@ -843,6 +1295,7 @@
close(15)
print *
+
end program meshfem2D
@@ -850,18 +1303,57 @@
! meshing subroutines
! *******************
+
!--- global node number
integer function num(i,j,nx)
- implicit none
+ implicit none
+
+ integer i,j,nx
+
+ num = j*(nx+1) + i + 1
+
+ end function num
+
- integer i,j,nx
+ !--- global node number (when ngnod==4).
+ integer function num_4(i,j,nx)
- num = j*(nx+1) + i + 1
+ implicit none
+
+ integer i,j,nx
+
+ num_4 = j*(nx+1) + i + 1
+
+ end function num_4
- end function num
+
+ !--- global node number (when ngnod==9).
+ integer function num_9(i,j,nx,nz)
+
+ implicit none
+
+ integer i,j,nx,nz
+
+
+ if ( (mod(i,2) == 0) .and. (mod(j,2) == 0) ) then
+ num_9 = j/2 * (nx+1) + i/2 + 1
+ else
+ if ( mod(j,2) == 0 ) then
+ num_9 = (nx+1)*(nz+1) + j/2 * nx + ceiling(real(i)/real(2))
+ else
+ num_9 = (nx+1)*(nz+1) + nx*(nz+1) + floor(real(j)/real(2))*(nx*2+1) + i + 1
+
+ end if
+ end if
+
+
+ end function num_9
+
+
+
!--- spline to describe the interfaces
double precision function value_spline(x,xinterface,zinterface,coefs_interface,npoints_interface)
Added: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -0,0 +1,1559 @@
+module part_unstruct
+
+ implicit none
+
+ include './constants_unstruct.h'
+
+
+contains
+
+
+ !-----------------------------------------------
+ ! Read the mesh and storing it in array 'elmnts' (which is allocated here).
+ ! 'num_start' is used to have the numbering of the nodes starting at '0'.
+ ! 'nelmnts' is the number of elements, 'nnodes' is the number of nodes in the mesh.
+ !-----------------------------------------------
+ subroutine read_mesh(filename, nelmnts, elmnts, nnodes, num_start)
+
+ character(len=256), intent(in) :: filename
+ integer, intent(out) :: nelmnts
+ integer, intent(out) :: nnodes
+ integer, dimension(:), pointer :: elmnts
+ integer, intent(out) :: num_start
+
+ integer :: i
+
+ print *, trim(filename)
+
+ open(unit=990, file=trim(filename), form='formatted' , status='old', action='read')
+ read(990,*) nelmnts
+ allocate(elmnts(0:ESIZE*nelmnts-1))
+ do i = 0, nelmnts-1
+ read(990,*) elmnts(i*ESIZE), elmnts(i*ESIZE+1), elmnts(i*ESIZE+2), elmnts(i*ESIZE+3)
+
+ end do
+ close(990)
+
+ num_start = minval(elmnts)
+ elmnts(:) = elmnts(:) - num_start
+ nnodes = maxval(elmnts) + 1
+
+
+ end subroutine read_mesh
+
+
+
+ subroutine read_nodes_coords(filename, nnodes, nodes_coords)
+
+ character(len=256), intent(in) :: filename
+ integer, intent(out) :: nnodes
+ double precision, dimension(:,:), pointer :: nodes_coords
+
+ integer :: i
+
+ print *, trim(filename)
+
+ open(unit=991, file=trim(filename), form='formatted' , status='old', action='read')
+ read(991,*) nnodes
+ allocate(nodes_coords(2,nnodes))
+ do i = 1, nnodes
+ read(991,*) nodes_coords(1,i), nodes_coords(2,i)
+
+ end do
+ close(991)
+
+ end subroutine read_nodes_coords
+
+
+
+ subroutine read_mat(filename, nelmnts, num_material)
+
+ character(len=256), intent(in) :: filename
+ integer, intent(in) :: nelmnts
+ integer, dimension(1:nelmnts), intent(out) :: num_material
+
+ integer :: i
+
+ print *, trim(filename)
+
+ open(unit=992, file=trim(filename), form='formatted' , status='old', action='read')
+ do i = 1, nelmnts
+ read(992,*) num_material(i)
+
+ end do
+ close(992)
+
+ end subroutine read_mat
+
+
+
+ subroutine read_acoustic_surface(filename, nelem_acoustic_surface, acoustic_surface, &
+ nelmnts, num_material, ANISOTROPIC_MATERIAL, nb_materials, icodemat, cs, num_start)
+
+ include './constants.h'
+
+ character(len=256), intent(in) :: filename
+ integer, intent(out) :: nelem_acoustic_surface
+ integer, dimension(:,:), pointer :: acoustic_surface
+ integer, intent(in) :: nelmnts
+ integer, dimension(0:nelmnts-1) :: num_material
+ integer, intent(in) :: ANISOTROPIC_MATERIAL
+ integer, intent(in) :: nb_materials
+ integer, dimension(1:nb_materials), intent(in) :: icodemat
+ double precision, dimension(1:nb_materials), intent(in) :: cs
+ integer, intent(in) :: num_start
+
+
+ integer, dimension(:,:), allocatable :: acoustic_surface_tmp
+ integer :: nelmnts_surface
+ integer :: i
+ integer :: imaterial_number
+
+
+ open(unit=993, file=trim(filename), form='formatted' , status='old', action='read')
+ read(993,*) nelmnts_surface
+ print *, 'POY', nelmnts_surface
+
+ allocate(acoustic_surface_tmp(4,nelmnts_surface))
+
+ do i = 1, nelmnts_surface
+ read(993,*) acoustic_surface_tmp(1,i), acoustic_surface_tmp(2,i), acoustic_surface_tmp(3,i), acoustic_surface_tmp(4,i)
+
+ end do
+
+ close(993)
+ acoustic_surface_tmp(1,:) = acoustic_surface_tmp(1,:) - num_start
+ acoustic_surface_tmp(3,:) = acoustic_surface_tmp(3,:) - num_start
+ acoustic_surface_tmp(4,:) = acoustic_surface_tmp(4,:) - num_start
+
+ nelem_acoustic_surface = 0
+ do i = 1, nelmnts_surface
+ imaterial_number = num_material(acoustic_surface_tmp(1,i))
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+
+ end if
+ end do
+
+ allocate(acoustic_surface(4,nelem_acoustic_surface))
+
+ nelem_acoustic_surface = 0
+ do i = 1, nelmnts_surface
+ imaterial_number = num_material(acoustic_surface_tmp(1,i))
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ acoustic_surface(:,nelem_acoustic_surface) = acoustic_surface_tmp(:,i)
+ end if
+ end do
+
+
+ end subroutine read_acoustic_surface
+
+
+
+ subroutine read_abs_surface(filename, nelemabs, abs_surface, num_start)
+
+ include './constants.h'
+
+ character(len=256), intent(in) :: filename
+ integer, intent(out) :: nelemabs
+ integer, dimension(:,:), pointer :: abs_surface
+ integer, intent(in) :: num_start
+
+
+ integer :: i
+
+
+ open(unit=994, file=trim(filename), form='formatted' , status='old', action='read')
+ read(994,*) nelemabs
+
+ allocate(abs_surface(4,nelemabs))
+
+ do i = 1, nelemabs
+ read(994,*) abs_surface(1,i), abs_surface(2,i), abs_surface(3,i), abs_surface(4,i)
+
+ end do
+
+ close(994)
+
+ abs_surface(1,:) = abs_surface(1,:) - num_start
+ abs_surface(3,:) = abs_surface(3,:) - num_start
+ abs_surface(4,:) = abs_surface(4,:) - num_start
+
+
+ end subroutine read_abs_surface
+
+
+
+ !*************************************************************************
+ ! creation du graphe dual (adjacence des elements definie par au moins un noeud commun )
+ !**************************************************************************
+ ! Version qui construit les adjacences en commencant par dresser la liste des elements pour chaque noeud (vertices).
+ ! On passe en plus en argument le tableau alloue qui contiendra la liste des elements par noeud, accompagne de la liste du
+ ! nombre d'elements qui contiennent chaque noeud (c'est pour pouvoir faire l'allocation avant sans connaitre la taille exacte).
+ subroutine mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts, ncommonnodes)
+
+ integer, intent(in) :: nelmnts
+ integer, intent(in) :: nnodes
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(:), pointer :: xadj
+ integer, dimension(:), pointer :: adjncy
+ integer, dimension(:), pointer :: nnodes_elmnts
+ integer, dimension(:), pointer :: nodes_elmnts
+ integer, intent(in) :: ncommonnodes
+
+ integer :: i, j, k, l, m, nb_edges
+ logical :: is_neighbour
+ integer :: num_node, n
+ integer :: elem_base, elem_target
+ integer :: connectivity
+
+
+ allocate(xadj(0:nelmnts))
+ xadj(:) = 0
+ allocate(adjncy(0:max_neighbour*nelmnts-1))
+ adjncy(:) = 0
+ allocate(nnodes_elmnts(0:nnodes-1))
+ nnodes_elmnts(:) = 0
+ allocate(nodes_elmnts(0:nsize*nnodes-1))
+ nodes_elmnts(:) = 0
+
+ nb_edges = 0
+
+ ! remplissage de la liste des elements par noeuds
+ do i = 0, esize*nelmnts-1
+ nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/esize
+ nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
+
+ end do
+
+ print *, 'nnodes_elmnts'
+ !print *, nnodes_elmnts
+
+ ! pour chaque noeud, remplissage du tableau des adjacences,
+ ! avec verification qu'un element ne soit pas compte plusieurs fois comme voisins d'un autre.
+ do j = 0, nnodes-1
+ do k = 0, nnodes_elmnts(j)-1
+ do l = k+1, nnodes_elmnts(j)-1
+
+ connectivity = 0
+ elem_base = nodes_elmnts(k+j*nsize)
+ elem_target = nodes_elmnts(l+j*nsize)
+ do n = 1, esize
+ num_node = elmnts(esize*elem_base+n-1)
+ do m = 0, nnodes_elmnts(num_node)-1
+ if ( nodes_elmnts(m+num_node*nsize) == elem_target ) then
+ connectivity = connectivity + 1
+ end if
+ end do
+ end do
+
+ if ( connectivity >= ncommonnodes) then
+
+ is_neighbour = .false.
+
+ do m = 0, xadj(nodes_elmnts(k+j*nsize))
+ if ( .not.is_neighbour ) then
+ if ( adjncy(nodes_elmnts(k+j*nsize)*max_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)*max_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
+ adjncy(nodes_elmnts(l+j*nsize)*max_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
+ end if
+ end if
+ end do
+ end do
+ end do
+
+ ! maintenant que la liste d'adjacence est contruite, on va refondre les tableaux dadjncy et dxadj pour qu'ils
+ ! soient au bon format
+ do i = 0, nelmnts-1
+ k = xadj(i)
+ xadj(i) = nb_edges
+ do j = 0, k-1
+ adjncy(nb_edges) = adjncy(i*max_neighbour+j)
+ nb_edges = nb_edges + 1
+ end do
+ end do
+
+ xadj(nelmnts) = nb_edges
+
+
+ end subroutine mesh2dual_ncommonnodes
+
+
+
+ subroutine read_weights(nelmnts, vwgt, nb_edges, adjwgt)
+
+ integer, intent(in) :: nelmnts, nb_edges
+ integer, dimension(:), pointer :: vwgt, adjwgt
+
+ allocate(vwgt(0:nelmnts-1))
+ allocate(adjwgt(0:nb_edges-1))
+
+ vwgt(:) = 1
+ adjwgt(:) = 1
+
+
+ end subroutine read_weights
+
+
+
+ !*************************************************************************
+ ! construction du tableau de correspondance entre numerotation globale et locale (pour les elements)
+ !*************************************************************************/
+ subroutine Construct_glob2loc_elmnts(nelmnts, part, nparts, glob2loc_elmnts)
+
+ integer, intent(in) :: nelmnts, nparts
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(:), pointer :: glob2loc_elmnts
+
+ integer :: num_glob, num_part
+ integer, dimension(0:nparts-1) :: num_loc
+
+
+ allocate(glob2loc_elmnts(0:nelmnts-1))
+
+ do num_part = 0, nparts-1
+ num_loc(num_part) = 0
+
+ end do
+
+ do num_glob = 0, nelmnts-1
+ num_part = part(num_glob)
+ glob2loc_elmnts(num_glob) = num_loc(num_part)
+ num_loc(num_part) = num_loc(num_part) + 1
+
+ end do
+
+
+ end subroutine Construct_glob2loc_elmnts
+
+
+
+ !*************************************************************************
+ ! construction du tableau de correspondance entre numerotation globale et locale (pour les noeuds)
+ !*************************************************************************/
+ subroutine Construct_glob2loc_nodes(nelmnts, nnodes, nnodes_elmnts, nodes_elmnts, part, nparts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
+
+ integer, intent(in) :: nelmnts, nnodes, nparts
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(0:nnodes-1), intent(in) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1), intent(in) :: nodes_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer :: num_node
+ integer :: el
+ integer :: num_part
+ integer :: size_glob2loc_nodes
+ integer, dimension(0:nparts-1) :: parts_node
+ integer, dimension(0:nparts-1) :: num_parts
+
+
+ allocate(glob2loc_nodes_nparts(0:nnodes))
+
+ size_glob2loc_nodes = 0
+
+!!$ do num_part = 0, nparts-1
+!!$ parts_node(num_part) = 0
+!!$ end do
+ parts_node(:) = 0
+
+
+ do num_node = 0, nnodes-1
+ 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
+
+ end do
+
+
+ glob2loc_nodes_nparts(nnodes) = size_glob2loc_nodes
+ !print *, glob2loc_nodes_nparts
+ !print *, glob2loc_nodes_nparts(nnodes)
+
+
+ allocate(glob2loc_nodes_parts(0:glob2loc_nodes_nparts(nnodes)-1))
+ allocate(glob2loc_nodes(0:glob2loc_nodes_nparts(nnodes)-1))
+
+ glob2loc_nodes(0) = 0
+
+!!$ do num_part = 0, nparts-1
+!!$ parts_node(num_part) = 0
+!!$ num_parts(num_part) = 0
+!!$
+!!$ end do
+ parts_node(:) = 0
+ num_parts(:) = 0
+ size_glob2loc_nodes = 0
+
+
+ do num_node = 0, nnodes-1
+ !print *, 'num_node', num_node
+ do el = 0, nnodes_elmnts(num_node)-1
+ !print *, 'el_', el
+ parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
+
+ end do
+ !print *, 'PPP'
+ do num_part = 0, nparts-1
+ !print *, 'num_part', num_part
+ !print *, 'parts_nodes(num_part)', parts_node(num_part)
+
+ if ( parts_node(num_part) == 1 ) then
+ glob2loc_nodes_parts(size_glob2loc_nodes) = num_part
+ glob2loc_nodes(size_glob2loc_nodes) = num_parts(num_part)
+ size_glob2loc_nodes = size_glob2loc_nodes + 1
+ num_parts(num_part) = num_parts(num_part) + 1
+ parts_node(num_part) = 0
+ end if
+ !print *, 'num_part__', num_part
+
+ end do
+ end do
+
+
+ end subroutine Construct_glob2loc_nodes
+
+
+
+ subroutine Construct_interfaces(nelmnts, nparts, part, elmnts, xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, nb_materials, cs_material, num_material)
+
+ include 'constants.h'
+
+ integer, intent(in) :: nelmnts, nparts
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(0:nelmnts), intent(in) :: xadj
+ integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(:),pointer :: tab_size_interfaces, tab_interfaces
+ integer, intent(out) :: ninterfaces
+ integer, dimension(1:nelmnts), intent(in) :: num_material
+ double precision, dimension(1:nb_materials), intent(in) :: cs_material
+ integer, intent(in) :: nb_materials
+
+
+ integer :: num_part, num_part_bis, el, el_adj, num_interface, num_edge, ncommon_nodes, &
+ num_node, num_node_bis
+ integer :: i, j
+ logical :: is_acoustic_el, is_acoustic_el_adj
+
+ print *, '00'
+ ninterfaces = 0
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ ninterfaces = ninterfaces + 1
+ end do
+ end do
+
+ allocate(tab_size_interfaces(0:ninterfaces))
+ tab_size_interfaces(:) = 0
+
+ print *, 'num_interface', ninterfaces
+
+ num_interface = 0
+ num_edge = 0
+ print *, '01'
+
+ do num_part = 0, nparts-1
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ if ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
+ do el_adj = xadj(el), xadj(el+1)-1
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ 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) ) then
+ num_edge = num_edge + 1
+
+ end if
+ end do
+ end if
+ end do
+ print *, 'num_interface', num_interface
+ tab_size_interfaces(num_interface+1) = tab_size_interfaces(num_interface) + num_edge
+ num_edge = 0
+ num_interface = num_interface + 1
+
+ end do
+ end do
+ print *, 'POUM',tab_size_interfaces(num_interface)
+
+ num_interface = 0
+ num_edge = 0
+
+ print *, '02'
+
+ ! le 5 vient du fait que l'on represente l'interface comme un tableau avec 5 valeurs: element1-element2-valence-noeud1-noeud2
+ allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*5-1)))
+ tab_interfaces(:) = 0
+ print *, '03'
+
+ do num_part = 0, nparts-1
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ if ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
+ do el_adj = xadj(el), xadj(el+1)-1
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ 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) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+0) = el
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+1) = adjncy(el_adj)
+ ncommon_nodes = 0
+ do num_node = 0, 4-1
+ do num_node_bis = 0, 4-1
+ if ( elmnts(el*esize+num_node) == elmnts(adjncy(el_adj)*esize+num_node_bis) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+3+ncommon_nodes) &
+ = elmnts(el*esize+num_node)
+ ncommon_nodes = ncommon_nodes + 1
+ end if
+ end do
+ end do
+ if ( ncommon_nodes > 0 ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+2) = ncommon_nodes
+ else
+ print *, "Erreur lors de la construction des interfaces!!!", ncommon_nodes
+ end if
+ num_edge = num_edge + 1
+ end if
+ end do
+ end if
+
+ end do
+ num_edge = 0
+ num_interface = num_interface + 1
+ end do
+ end do
+
+
+ end subroutine Construct_interfaces
+
+
+
+ subroutine write_glob2loc_nodes_database(IIN_database, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, num_phase)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: nnodes, iproc, num_phase
+ integer, intent(inout) :: npgeo
+
+ double precision, dimension(2,nnodes) :: nodes_coords
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer :: i, j
+
+ if ( num_phase == 1 ) then
+ npgeo = 0
+
+ do i = 0, nnodes-1
+ do j = glob2loc_nodes_nparts(i), glob2loc_nodes_nparts(i+1)-1
+ if ( glob2loc_nodes_parts(j) == iproc ) then
+ npgeo = npgeo + 1
+
+ end if
+ !write(992,*) i, glob2loc_nodes_parts(j), glob2loc_nodes(j)
+
+ end do
+ end do
+ else
+ 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)
+ end if
+ end do
+ end do
+ end if
+
+ end subroutine Write_glob2loc_nodes_database
+
+
+
+ subroutine write_partition_database(IIN_database, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, num_modele, ngnod, num_phase)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: nelmnts, num_phase, iproc
+ integer, intent(inout) :: nspec
+ integer, dimension(:), pointer :: part, elmnts, glob2loc_elmnts
+ integer, dimension(:) :: num_modele
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, intent(in) :: ngnod
+
+ integer :: i,j,k
+ integer, dimension(0:ngnod-1) :: loc_nodes
+
+
+ if ( num_phase == 1 ) then
+ nspec = 0
+
+ do i = 0, nelmnts-1
+ if ( part(i) == iproc ) then
+ nspec = nspec + 1
+
+ end if
+ end do
+
+ else
+ do i = 0, nelmnts-1
+ if ( part(i) == iproc ) then
+
+ do j = 0, ngnod-1
+ do k = glob2loc_nodes_nparts(elmnts(i*ngnod+j)), glob2loc_nodes_nparts(elmnts(i*ngnod+j)+1)-1
+
+ if ( glob2loc_nodes_parts(k) == iproc ) then
+ loc_nodes(j) = glob2loc_nodes(k)
+
+ end if
+ end do
+
+ end do
+ write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(i+1), (loc_nodes(k)+1, k=0,ngnod-1)
+ end if
+ end do
+ end if
+
+
+ end subroutine write_partition_database
+
+
+
+ subroutine Write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, nparts, iproc, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, num_phase)
+
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer, intent(in) :: nparts
+ integer, intent(in) :: ninterfaces
+ integer, intent(inout) :: my_ninterface
+ integer, dimension(:), pointer :: tab_size_interfaces
+ integer, dimension(:), pointer :: tab_interfaces
+ integer, dimension(0:ninterfaces-1), intent(inout) :: my_interfaces
+ integer, dimension(0:ninterfaces-1), intent(inout) :: my_nb_interfaces
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+
+ integer, dimension(2) :: local_nodes
+ integer :: local_elmnt
+ integer :: num_phase
+
+ integer :: i, j, k, l
+ integer :: num_interface
+
+
+ num_interface = 0
+
+ if ( num_phase == 1 ) then
+
+ my_interfaces(:) = 0
+ my_nb_interfaces(:) = 0
+
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ if ( (tab_size_interfaces(num_interface) < tab_size_interfaces(num_interface+1)) .and. &
+ (i == iproc .or. j == iproc) ) then
+ my_interfaces(num_interface) = 1
+ 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
+ end do
+ my_ninterface = sum(my_interfaces(:))
+
+ else
+
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ if ( my_interfaces(num_interface) == 1 ) then
+ if ( i == iproc ) then
+ write(IIN_database,*) j, my_nb_interfaces(num_interface)
+ else
+ write(IIN_database,*) i, my_nb_interfaces(num_interface)
+ end if
+
+ do k = tab_size_interfaces(num_interface), tab_size_interfaces(num_interface+1)-1
+ if ( i == iproc ) then
+ local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+0))+1
+ else
+ local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+1))+1
+ end if
+
+ if ( tab_interfaces(k*5+2) == 1 ) then
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
+ glob2loc_nodes_nparts(tab_interfaces(k*5+3)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(1) = glob2loc_nodes(l)+1
+ end if
+ end do
+
+ write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), local_nodes(1), -1
+ else
+ if ( tab_interfaces(k*5+2) == 2 ) then
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
+ glob2loc_nodes_nparts(tab_interfaces(k*5+3)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(1) = glob2loc_nodes(l)+1
+ end if
+ end do
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*5+4)), &
+ glob2loc_nodes_nparts(tab_interfaces(k*5+4)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(2) = glob2loc_nodes(l)+1
+ end if
+ end do
+ write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), local_nodes(1), local_nodes(2)
+ else
+ write(IIN_database,*) "erreur_write_interface_", tab_interfaces(k*5+2)
+ end if
+ end if
+ end do
+
+ end if
+
+ num_interface = num_interface + 1
+ end do
+ end do
+
+ end if
+
+
+ end subroutine Write_interfaces_database
+
+
+
+ subroutine Write_surface_database(IIN_database, nsurface, surface, &
+ nsurface_loc, nparts, iproc, glob2loc_elmnts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, num_phase)
+
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer, intent(in) :: nparts
+ integer :: nsurface
+ integer :: nsurface_loc
+ integer, dimension(:,:), pointer :: surface
+
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, dimension(:), pointer :: part
+
+
+ integer, dimension(2) :: local_nodes
+ integer :: local_elmnt
+ integer :: num_phase
+
+ integer :: i, j, k, l
+ integer :: num_surface
+
+
+
+ if ( num_phase == 1 ) then
+
+ nsurface_loc = 0
+
+ do i = 1, nsurface
+ if ( part(surface(1,i)) == iproc ) then
+ nsurface_loc = nsurface_loc + 1
+
+ end if
+ end do
+
+ else
+
+ nsurface_loc = 0
+
+ do i = 1, nsurface
+ if ( part(surface(1,i)) == iproc ) then
+ nsurface_loc = nsurface_loc + 1
+
+ local_elmnt = glob2loc_elmnts(surface(1,i)) + 1
+
+ if ( surface(2,i) == 1 ) then
+ do l = glob2loc_nodes_nparts(surface(3,i)), &
+ glob2loc_nodes_nparts(surface(3,i)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(1) = glob2loc_nodes(l)+1
+ end if
+ end do
+
+ write(IIN_database,*) local_elmnt, surface(2,i), local_nodes(1), -1
+ end if
+ if ( surface(2,i) == 2 ) then
+ do l = glob2loc_nodes_nparts(surface(3,i)), &
+ glob2loc_nodes_nparts(surface(3,i)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(1) = glob2loc_nodes(l)+1
+ end if
+ end do
+ do l = glob2loc_nodes_nparts(surface(4,i)), &
+ glob2loc_nodes_nparts(surface(4,i)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(2) = glob2loc_nodes(l)+1
+ end if
+ end do
+
+ write(IIN_database,*) local_elmnt, surface(2,i), local_nodes(1), local_nodes(2)
+ end if
+
+ end if
+
+ end do
+
+ end if
+
+
+ end subroutine Write_surface_database
+
+
+
+ subroutine merge_abs_boundaries(nelemabs, nelemabs_merge, abs_surface, abs_surface_char, abs_surface_merge, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ nodes_coords, &
+ nedges_coupled, edges_coupled, nb_materials, cs_material, num_material, &
+ nelmnts, &
+ elmnts, ngnod)
+
+ implicit none
+ include 'constants.h'
+
+
+ integer, intent(inout) :: nelemabs
+ integer, intent(out) :: nelemabs_merge
+ integer, dimension(:,:), pointer :: abs_surface
+ logical, dimension(:,:), pointer :: abs_surface_char
+ integer, dimension(:), pointer :: abs_surface_merge
+ integer, dimension(:), pointer :: elmnts
+ integer, intent(in) :: ngnod
+ integer, dimension(:), pointer :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right
+ double precision, dimension(:,:), pointer :: nodes_coords
+ integer :: nedges_coupled
+ integer, dimension(:,:), pointer :: edges_coupled
+ integer :: nb_materials
+ double precision, dimension(nb_materials), intent(in) :: cs_material
+ integer, dimension(1:nelmnts), intent(in) :: num_material
+ integer :: nelmnts
+
+
+ logical, dimension(nb_materials) :: is_acoustic
+ integer :: num_edge, nedge_bound
+ integer :: num_edge_bis
+ integer :: match
+ integer :: nb_elmnts_abs
+ integer :: i
+ integer :: temp
+ integer :: common_node, other_node
+ integer :: modified_edge_elmnt
+ integer :: modified_edge
+ double precision :: vect_product
+ integer :: iedge, inode1, inode2
+
+
+
+ allocate(abs_surface_char(4,nelemabs))
+ allocate(abs_surface_merge(nelemabs))
+ abs_surface_char(:,:) = .false.
+ abs_surface_merge(:) = -1
+
+ nedge_bound = nelemabs
+ nb_elmnts_abs = 0
+
+ do num_edge = 1, nedge_bound
+
+ match = 0
+ do i = 1, nb_elmnts_abs
+ if ( abs_surface(1,num_edge) == abs_surface_merge(i) ) then
+ match = i
+ exit
+ end if
+ end do
+
+ if ( match == 0 ) then
+ nb_elmnts_abs = nb_elmnts_abs + 1
+ match = nb_elmnts_abs
+ end if
+ !print *, 'match', num_edge,match
+
+ abs_surface_merge(match) = abs_surface(1,num_edge)
+
+
+ if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1)) ) then
+ abs_surface_char(1,match) = .true.
+
+ end if
+
+ if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
+ abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1)) ) then
+ temp = abs_surface(4,num_edge)
+ abs_surface(4,num_edge) = abs_surface(3,num_edge)
+ abs_surface(3,num_edge) = temp
+ abs_surface_char(1,match) = .true.
+
+ end if
+
+ if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
+ abs_surface_char(4,match) = .true.
+
+ end if
+
+ if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
+ abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
+ temp = abs_surface(4,num_edge)
+ abs_surface(4,num_edge) = abs_surface(3,num_edge)
+ abs_surface(3,num_edge) = temp
+ abs_surface_char(4,match) = .true.
+
+ end if
+
+ if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2)) ) then
+ abs_surface_char(2,match) = .true.
+
+ end if
+
+ if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1) .and. &
+ abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2)) ) then
+ temp = abs_surface(4,num_edge)
+ abs_surface(4,num_edge) = abs_surface(3,num_edge)
+ abs_surface(3,num_edge) = temp
+ abs_surface_char(2,match) = .true.
+
+ end if
+
+ if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
+ temp = abs_surface(4,num_edge)
+ abs_surface(4,num_edge) = abs_surface(3,num_edge)
+ abs_surface(3,num_edge) = temp
+ abs_surface_char(3,match) = .true.
+
+ end if
+
+ if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2) .and. &
+ abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
+ abs_surface_char(3,match) = .true.
+
+ end if
+
+ end do
+
+ nelemabs_merge = nb_elmnts_abs
+
+
+ allocate(ibegin_bottom(nelemabs_merge))
+ allocate(iend_bottom(nelemabs_merge))
+ allocate(jbegin_right(nelemabs_merge))
+ allocate(jend_right(nelemabs_merge))
+ allocate(ibegin_top(nelemabs_merge))
+ allocate(iend_top(nelemabs_merge))
+ allocate(jbegin_left(nelemabs_merge))
+ allocate(jend_left(nelemabs_merge))
+
+ ibegin_bottom(:) = 1
+ jbegin_right(:) = 1
+ ibegin_top(:) = 1
+ jbegin_left(:) = 1
+ iend_bottom(:) = NGLLX
+ jend_right(:) = NGLLZ
+ iend_top(:) = NGLLX
+ jend_left(:) = NGLLZ
+
+
+!!$ do num_edge = 1, nedge_bound
+!!$ do num_edge_bis = 1, nedge_bound
+!!$
+!!$ modified_edge = 0
+!!$
+!!$ if ( abs_surface(1,num_edge) /= abs_surface(1,num_edge_bis) ) then
+!!$
+!!$ vect_product = ( (nodes_coords(1,abs_surface(4,num_edge)+1) - nodes_coords(1,abs_surface(3,num_edge)+1) ) &
+!!$ * (nodes_coords(2,abs_surface(4,num_edge_bis)+1) - nodes_coords(2,abs_surface(3,num_edge_bis)+1) ) &
+!!$ - (nodes_coords(2,abs_surface(4,num_edge)+1) - nodes_coords(2,abs_surface(3,num_edge)+1) ) &
+!!$ * (nodes_coords(1,abs_surface(4,num_edge_bis)+1) - nodes_coords(1,abs_surface(3,num_edge_bis)+1) ) )
+!!$ if ( abs(vect_product) > TINYVAL*1000000 ) then
+!!$ !if ( abs(vect_product) > 10 ) then
+!!$
+!!$ if ( abs_surface(3,num_edge) == abs_surface(3,num_edge_bis) ) then
+!!$ common_node = abs_surface(3,num_edge)
+!!$ if ( abs_surface(1,num_edge) < abs_surface(1,num_edge_bis) ) then
+!!$ modified_edge_elmnt = abs_surface(1,num_edge_bis)
+!!$ other_node = abs_surface(4,num_edge_bis)
+!!$ modified_edge = num_edge_bis
+!!$ else
+!!$ modified_edge_elmnt = abs_surface(1,num_edge)
+!!$ other_node = abs_surface(4,num_edge)
+!!$ modified_edge = num_edge
+!!$ end if
+!!$ end if
+!!$
+!!$ if ( abs_surface(3,num_edge) == abs_surface(4,num_edge_bis) ) then
+!!$ common_node = abs_surface(3,num_edge)
+!!$ if ( abs_surface(1,num_edge) < abs_surface(1,num_edge_bis) ) then
+!!$ modified_edge_elmnt = abs_surface(1,num_edge_bis)
+!!$ other_node = abs_surface(3,num_edge_bis)
+!!$ modified_edge = num_edge_bis
+!!$ else
+!!$ modified_edge_elmnt = abs_surface(1,num_edge)
+!!$ other_node = abs_surface(4,num_edge)
+!!$ modified_edge = num_edge
+!!$ end if
+!!$ end if
+!!$
+!!$ if ( abs_surface(4,num_edge) == abs_surface(3,num_edge_bis) ) then
+!!$ common_node = abs_surface(4,num_edge)
+!!$ if ( abs_surface(1,num_edge) < abs_surface(1,num_edge_bis) ) then
+!!$ modified_edge_elmnt = abs_surface(1,num_edge_bis)
+!!$ other_node = abs_surface(4,num_edge_bis)
+!!$ modified_edge = num_edge_bis
+!!$ else
+!!$ modified_edge_elmnt = abs_surface(1,num_edge)
+!!$ other_node = abs_surface(3,num_edge)
+!!$ modified_edge = num_edge
+!!$ end if
+!!$ end if
+!!$
+!!$ if ( abs_surface(4,num_edge) == abs_surface(4,num_edge_bis) ) then
+!!$ common_node = abs_surface(4,num_edge)
+!!$ if ( abs_surface(1,num_edge) < abs_surface(1,num_edge_bis) ) then
+!!$ modified_edge_elmnt = abs_surface(1,num_edge_bis)
+!!$ other_node = abs_surface(3,num_edge_bis)
+!!$ modified_edge = num_edge_bis
+!!$ else
+!!$ modified_edge_elmnt = abs_surface(1,num_edge)
+!!$ other_node = abs_surface(3,num_edge)
+!!$ modified_edge = num_edge
+!!$ end if
+!!$
+!!$ end if
+!!$
+!!$
+!!$ if ( modified_edge /= 0 ) then
+!!$ print *, 'SSSSSSS', nodes_coords(1,common_node+1), nodes_coords(2,common_node+1), &
+!!$ common_node, modified_edge, modified_edge_elmnt, other_node
+!!$ match = 0
+!!$ do i = 1, nelemabs_merge
+!!$ if ( abs_surface(1,modified_edge) == abs_surface_merge(i) ) then
+!!$ match = i
+!!$ exit
+!!$ end if
+!!$ end do
+!!$
+!!$ if ( abs_surface(3,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
+!!$ abs_surface(4,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+1) ) then
+!!$ if ( common_node == abs_surface(3,modified_edge) ) then
+!!$ ibegin_bottom(match) = 2
+!!$ else
+!!$ iend_bottom(match) = NGLLX - 1
+!!$ end if
+!!$ end if
+!!$
+!!$ if ( abs_surface(3,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+1) .and. &
+!!$ abs_surface(4,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+2) ) then
+!!$ if ( common_node == abs_surface(3,modified_edge) ) then
+!!$ jbegin_right(match) = 2
+!!$ else
+!!$ jend_right(match) = NGLLZ - 1
+!!$ end if
+!!$ end if
+!!$
+!!$ if ( abs_surface(3,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+3) .and. &
+!!$ abs_surface(4,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+2) ) then
+!!$ if ( common_node == abs_surface(3,modified_edge) ) then
+!!$ ibegin_top(match) = 2
+!!$ else
+!!$ iend_top(match) = NGLLX - 1
+!!$ end if
+!!$ end if
+!!$
+!!$ if ( abs_surface(3,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
+!!$ abs_surface(4,modified_edge) == elmnts(ngnod*abs_surface_merge(match)+3) ) then
+!!$ if ( common_node == abs_surface(3,modified_edge) ) then
+!!$ jbegin_left(match) = 2
+!!$ else
+!!$ jend_left(match) = NGLLZ - 1
+!!$ end if
+!!$ end if
+!!$
+!!$ end if
+!!$
+!!$ end if
+!!$
+!!$ end if
+!!$
+!!$ end do
+!!$
+!!$ end do
+
+
+ is_acoustic(:) = .false.
+ do i = 1, nb_materials
+ if (cs_material(i) < TINYVAL) then
+ is_acoustic(i) = .true.
+ end if
+
+ end do
+
+
+ do num_edge = 1, nedge_bound
+
+ match = 0
+ do i = 1, nelemabs_merge
+ if ( abs_surface(1,num_edge) == abs_surface_merge(i) ) then
+ match = i
+ exit
+ end if
+ end do
+
+ if ( is_acoustic(num_material(abs_surface(1,num_edge)+1)) ) then
+
+ do iedge = 1, nedges_coupled
+
+ do inode1 = 0, 3
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*edges_coupled(1,iedge)+inode1) ) then
+ do inode2 = 0, 3
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*edges_coupled(2,iedge)+inode2) ) then
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) ) then
+ ibegin_bottom(match) = 2
+
+ end if
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ jbegin_right(match) = 2
+
+ end if
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ ibegin_top(match) = 2
+
+ end if
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) ) then
+ jbegin_left(match) = 2
+
+ end if
+
+ end if
+ end do
+
+ end if
+
+ if ( abs_surface(4,num_edge) == elmnts(ngnod*edges_coupled(1,iedge)+inode1) ) then
+ do inode2 = 0, 3
+ if ( abs_surface(4,num_edge) == elmnts(ngnod*edges_coupled(2,iedge)+inode2) ) then
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) ) then
+ iend_bottom(match) = NGLLX - 1
+
+ end if
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ jend_right(match) = NGLLZ - 1
+
+ end if
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ iend_top(match) = NGLLX - 1
+
+ end if
+ if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) ) then
+ jend_left(match) = NGLLZ - 1
+
+ end if
+ end if
+ end do
+
+ end if
+
+ end do
+
+
+ end do
+
+ end if
+
+ end do
+
+!!$print *, 'AAAAAAAAa'
+!!$ do iedge = 1, nedges_coupled
+!!$ print *, nodes_coords(1,elmnts(ngnod*edges_coupled(1,iedge))+1), &
+!!$ nodes_coords(2,elmnts(ngnod*edges_coupled(1,iedge))+1), 0
+!!$ print *, nodes_coords(1,elmnts(ngnod*edges_coupled(2,iedge))+1), &
+!!$ nodes_coords(2,elmnts(ngnod*edges_coupled(2,iedge))+1), 0
+!!$
+!!$ end do
+!!$print *, 'AAAAAAAAa'
+
+ end subroutine merge_abs_boundaries
+
+
+
+ subroutine write_abs_merge_database(IIN_database, nelemabs_merge, nelemabs_loc, &
+ abs_surface_char, abs_surface_merge, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ glob2loc_elmnts, part, iproc, num_phase)
+
+ implicit none
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: nelemabs_merge
+ integer, intent(inout) :: nelemabs_loc
+ logical, dimension(:,:), pointer :: abs_surface_char
+ integer, dimension(:), pointer :: abs_surface_merge
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: part
+ integer, intent(in) :: iproc
+ integer, intent(in) :: num_phase
+ integer, dimension(:), pointer :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right
+
+ integer :: i
+
+
+ if ( num_phase == 1 ) then
+ nelemabs_loc = 0
+ do i = 1, nelemabs_merge
+ if ( part(abs_surface_merge(i)) == iproc ) then
+ nelemabs_loc = nelemabs_loc + 1
+ end if
+ end do
+ else
+ do i = 1, nelemabs_merge
+ if ( part(abs_surface_merge(i)) == iproc ) then
+
+ write(IIN_database,*) glob2loc_elmnts(abs_surface_merge(i))+1, abs_surface_char(1,i), &
+ abs_surface_char(2,i), abs_surface_char(3,i), abs_surface_char(4,i), &
+ ibegin_bottom(i), iend_bottom(i), &
+ jbegin_right(i), jend_right(i), &
+ ibegin_top(i), iend_top(i), &
+ jbegin_left(i), jend_left(i)
+
+ end if
+
+ end do
+ end if
+
+
+ end subroutine write_abs_merge_database
+
+
+
+#ifdef USE_METIS
+ subroutine Part_metis(nelmnts, xadj, adjncy, vwgt, adjwgt, nparts, nb_edges, edgecut, part, metis_options)
+
+ integer, intent(in) :: nelmnts, nparts, nb_edges
+ integer, intent(inout) :: edgecut
+ integer, dimension(0:nelmnts), intent(in) :: xadj
+ integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(0:nelmnts-1), intent(in) :: vwgt
+ integer, dimension(0:nb_edges-1), intent(in) :: adjwgt
+ integer, dimension(:), pointer :: part
+ integer, dimension(0:4) :: metis_options
+
+ integer :: wgtflag
+ integer :: num_start
+
+ num_start = 0
+ wgtflag = 0
+
+ print *, 'avant', edgecut
+ call METIS_PartGraphRecursive(nelmnts, xadj(0), adjncy(0), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
+ metis_options, edgecut, part(0));
+ !call METIS_PartGraphVKway(nelmnts, xadj(0), adjncy(0), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
+ ! options, edgecut, part(0));
+ print *, 'apres', edgecut
+
+
+ end subroutine Part_metis
+#endif
+
+
+
+#ifdef USE_SCOTCH
+ subroutine Part_scotch(nelmnts, xadj, adjncy, vwgt, adjwgt, nparts, nedges, edgecut, part, scotch_strategy)
+
+ include 'scotchf.h'
+
+ integer, intent(in) :: nelmnts, nparts, nedges
+ integer, intent(inout) :: edgecut
+ integer, dimension(0:nelmnts), intent(in) :: xadj
+ integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(0:nelmnts-1), intent(in) :: vwgt
+ integer, dimension(:), pointer :: adjwgt
+ integer, dimension(:), pointer :: part
+
+ double precision, dimension(SCOTCH_GRAPHDIM) :: SCOTCHGRAPH
+ double precision, dimension(SCOTCH_STRATDIM) :: SCOTCHSTRAT
+ character(len=256), intent(in) :: scotch_strategy
+ integer :: IERR
+
+ integer :: wgtflag
+ integer, dimension(0:4) :: options
+ integer :: num_start
+
+
+ call scotchfstratinit (SCOTCHSTRAT(1), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot initialize strat'
+ STOP
+ ENDIF
+
+ !scotch_strategy="b{strat=x}"
+ !scotch_strategy='b{strat=h{pass=10000}}'
+ !scotch_strategy='b{strat=m{asc=b{strat=g{pass=100}},low=g{pass=100},type=h,vert=10}}'
+ !scotch_strategy='b{strat=f{bal=0.0001,move=1000000,pass=-1}}'
+ !scotch_strategy='b{strat=g}'
+ !scotch_strategy='g'
+ !scotch_strategy='b{strat=hf{move=1000}}'
+ !scotch_strategy='b{strat=m{vert=4,asc=h|g,low=f}}'
+ !scotch_strategy='b{strat=m{asc=g,low=g,rat=1.0,type=h,vert=4}}'
+ !scotch_strategy='b{strat=m{asc=g,low=f{bal=0.001,move=10000,pass=500},rat=0.5,type=h,vert=4}}'
+ call scotchfstratgraphmap (SCOTCHSTRAT(1), trim(scotch_strategy), IERR)
+ !call scotchfstratgraphbipart (SCOTCHSTRAT(1), trim(scotch_strategy), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot build strat'
+ STOP
+ ENDIF
+
+ CALL SCOTCHFGRAPHINIT (SCOTCHGRAPH (1), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot initialize graph'
+ STOP
+ ENDIF
+
+ CALL SCOTCHFGRAPHBUILD (SCOTCHGRAPH (1), 0, nelmnts, &
+ xadj (0), xadj (0), &
+ xadj (0), xadj (0), &
+ nedges, &
+ ! adjncy (0), adjncy (0), IERR)
+ adjncy (0), adjwgt (0), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot build graph'
+ STOP
+ ENDIF
+
+ CALL SCOTCHFGRAPHCHECK (SCOTCHGRAPH (1), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Invalid check'
+ STOP
+ ENDIF
+
+ call scotchfgraphpart (SCOTCHGRAPH (1), nparts, SCOTCHSTRAT(1), part(0), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot part graph'
+ STOP
+ ENDIF
+
+ CALL SCOTCHFGRAPHEXIT (SCOTCHGRAPH (1), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot destroy graph'
+ STOP
+ ENDIF
+
+ call scotchfstratexit (SCOTCHSTRAT(1), IERR)
+ IF (IERR .NE. 0) THEN
+ PRINT *, 'ERROR : MAIN : Cannot destroy strat'
+ STOP
+ ENDIF
+
+
+ end subroutine Part_scotch
+#endif
+
+
+
+subroutine acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, nb_materials, cs_material, num_material, &
+ nproc, part, nedges_coupled, edges_coupled)
+
+ implicit none
+ include 'constants.h'
+
+ integer, intent(in) :: nelmnts, nnodes, nproc, nb_materials
+ double precision, dimension(nb_materials), intent(in) :: cs_material
+ integer, dimension(1:nelmnts), intent(in) :: num_material
+ integer, dimension(:), pointer :: elmnts
+ integer, dimension(:), pointer :: part
+ integer, intent(out) :: nedges_coupled
+ integer, dimension(:,:), pointer :: edges_coupled
+
+
+ logical, dimension(nb_materials) :: is_acoustic
+ integer, dimension(:), pointer :: xadj
+ integer, dimension(:), pointer :: adjncy
+ integer, dimension(:), pointer :: nodes_elmnts
+ integer, dimension(:), pointer :: nnodes_elmnts
+
+ integer :: i, num_edge
+ integer :: el, el_adj
+ logical :: is_repartitioned
+
+
+ integer :: aaa
+
+
+ is_acoustic(:) = .false.
+ do i = 1, nb_materials
+ if (cs_material(i) < TINYVAL) then
+ is_acoustic(i) = .true.
+ end if
+
+ end do
+
+
+ call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,2)
+
+!!$ aaa = 1600
+!!$ print *, 'POUM', xadj(aaa), xadj(aaa+1) - 1
+!!$ do i = xadj(aaa), xadj(aaa+1) -1
+!!$ print *, i, adjncy(i)
+!!$ end do
+!!$ stop
+
+ nedges_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_acoustic(num_material(el+1)) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if ( .not. is_acoustic(num_material(adjncy(el_adj)+1)) ) then
+ nedges_coupled = nedges_coupled + 1
+ end if
+
+ end do
+ end if
+ end do
+
+ print *, 'nedges_coupled', nedges_coupled
+
+ allocate(edges_coupled(2,nedges_coupled))
+
+ nedges_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_acoustic(num_material(el+1)) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if ( .not. is_acoustic(num_material(adjncy(el_adj)+1)) ) then
+ nedges_coupled = nedges_coupled + 1
+ edges_coupled(1,nedges_coupled) = el
+ edges_coupled(2,nedges_coupled) = adjncy(el_adj)
+ end if
+
+ end do
+ end if
+ end do
+
+ do i = 1, nedges_coupled*nproc
+ is_repartitioned = .false.
+ do num_edge = 1, nedges_coupled
+ if ( part(edges_coupled(1,num_edge)) /= part(edges_coupled(2,num_edge)) ) then
+ if ( part(edges_coupled(1,num_edge)) < part(edges_coupled(2,num_edge)) ) then
+ part(edges_coupled(2,num_edge)) = part(edges_coupled(1,num_edge))
+ else
+ part(edges_coupled(1,num_edge)) = part(edges_coupled(2,num_edge))
+ end if
+ is_repartitioned = .true.
+ end if
+
+ end do
+ if ( .not. is_repartitioned ) then
+ exit
+ end if
+ end do
+
+
+end subroutine acoustic_elastic_repartitioning
+
+
+
+subroutine write_acoustic_elastic_edges_database(IIN_database, nedges_coupled, nedges_coupled_loc, &
+ edges_coupled, glob2loc_elmnts, nelmnts,part, iproc, num_phase)
+
+ implicit none
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: nedges_coupled
+ integer, intent(inout) :: nedges_coupled_loc
+ integer, dimension(:,:), pointer :: edges_coupled
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, intent(in) :: nelmnts
+ integer, dimension(:), pointer :: part
+ integer, intent(in) :: iproc
+ integer, intent(in) :: num_phase
+
+
+ integer :: i
+
+
+ if ( num_phase == 1 ) then
+ nedges_coupled_loc = 0
+ do i = 1, nedges_coupled
+ if ( part(edges_coupled(1,i)) == iproc ) then
+ nedges_coupled_loc = nedges_coupled_loc + 1
+ end if
+ end do
+ else
+ do i = 1, nedges_coupled
+ if ( part(edges_coupled(1,i)) == iproc ) then
+ write(IIN_database,*) glob2loc_elmnts(edges_coupled(1,i))+1, glob2loc_elmnts(edges_coupled(2,i))+1
+
+ end if
+
+ end do
+ end if
+
+
+end subroutine write_acoustic_elastic_edges_database
+
+
+
+end module part_unstruct
Added: seismo/2D/SPECFEM2D/trunk/scotchf.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/scotchf.h 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/scotchf.h 2007-12-07 23:53:26 UTC (rev 8524)
@@ -0,0 +1,14 @@
+ INTEGER SCOTCH_ARCHDIM
+ INTEGER SCOTCH_GEOMDIM
+ INTEGER SCOTCH_GRAPHDIM
+ INTEGER SCOTCH_MAPDIM
+ INTEGER SCOTCH_MESHDIM
+ INTEGER SCOTCH_ORDERDIM
+ INTEGER SCOTCH_STRATDIM
+ PARAMETER (SCOTCH_ARCHDIM = 4)
+ PARAMETER (SCOTCH_GEOMDIM = 2)
+ PARAMETER (SCOTCH_GRAPHDIM = 12)
+ PARAMETER (SCOTCH_MAPDIM = 10)
+ PARAMETER (SCOTCH_MESHDIM = 15)
+ PARAMETER (SCOTCH_ORDERDIM = 11)
+ PARAMETER (SCOTCH_STRATDIM = 1)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -91,6 +91,9 @@
implicit none
include "constants.h"
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
character(len=80) datlin
@@ -112,11 +115,7 @@
! pressure in an element
double precision, dimension(NGLLX,NGLLX) :: pressure_element
-! to write seismograms in single precision SEP and double precision binary format
- real(kind=4), dimension(:), allocatable :: buffer_binary_single
- double precision, dimension(:), allocatable :: buffer_binary_double
-
- integer :: i,j,k,it,irec,ipoin,ip,id,nbpoin,inump,n,ispec,iedge,npoin,npgeo,iglob
+ integer :: i,j,k,l,it,irec,ipoin,ip,id,nbpoin,inump,n,ispec,iedge,npoin,npgeo,iglob
logical :: anyabs
double precision :: dxd,dzd,valux,valuz,hlagrange,rhol,cosrot,sinrot,xi,gamma,x,z
@@ -143,7 +142,6 @@
! for acoustic medium
double precision, dimension(:), allocatable :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
-
double precision, dimension(:), allocatable :: rmass_inverse_elastic,rmass_inverse_acoustic,density,displread,velocread,accelread
double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
@@ -158,7 +156,7 @@
integer, dimension(:), allocatable :: kmato,numabs,ispecnum_acoustic_surface,iedgenum_acoustic_surface, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top,jbegin_left,jend_left,jbegin_right,jend_right
- integer ispec_selected_source,iglob_source,ix_source,iz_source
+ integer ispec_selected_source,iglob_source,ix_source,iz_source,is_proc_source,nb_proc_source
double precision a,displnorm_all
double precision, dimension(:), allocatable :: source_time_function
double precision, external :: erf
@@ -207,6 +205,23 @@
integer, dimension(:,:), allocatable :: iglob_image_color,copy_iglob_image_color
double precision, dimension(:,:), allocatable :: image_color_data
+ double precision :: xmin_color_image_loc, xmax_color_image_loc, zmin_color_image_loc, &
+ zmax_color_image_loc
+ integer :: min_i, min_j, max_i, max_j
+ integer :: nb_pixel_loc
+ integer, dimension(:), allocatable :: nb_pixel_per_proc
+ double precision :: i_coord, j_coord
+ double precision, dimension(2,4) :: elmnt_coords
+ integer, dimension(:), allocatable :: num_pixel_loc
+ integer, dimension(:,:), allocatable :: num_pixel_recv
+ double precision, dimension(:), allocatable :: data_pixel_recv
+ double precision, dimension(:), allocatable :: data_pixel_send
+ logical :: pixel_is_in
+ double precision :: dist_pixel, dist_min_pixel
+#ifdef USE_MPI
+ integer, dimension(MPI_STATUS_SIZE) :: request_mpi_status
+#endif
+
! timing information for the stations
character(len=MAX_LENGTH_STATION_NAME), allocatable, dimension(:) :: station_name
character(len=MAX_LENGTH_NETWORK_NAME), allocatable, dimension(:) :: network_name
@@ -228,15 +243,70 @@
integer, dimension(8) :: time_values
integer ihours,iminutes,iseconds,int_tCPU
double precision :: time_start,time_end,tCPU
+
+! for MPI and partitionning
+ integer :: ier
+ integer :: nproc
+ integer :: myrank
+ integer :: iproc
+ character(len=256) :: prname
+ integer :: ninterface
+ integer :: max_interface_size
+ integer, dimension(:), allocatable :: my_neighbours
+ integer, dimension(:), allocatable :: my_nelmnts_neighbours
+ integer, dimension(:,:,:), allocatable :: my_interfaces
+ integer, dimension(:,:), allocatable :: ibool_interfaces_acoustic,ibool_interfaces_elastic
+ integer, dimension(:), allocatable :: nibool_interfaces_acoustic,nibool_interfaces_elastic
+
+ integer :: ninterface_acoustic, ninterface_elastic
+ integer, dimension(:), allocatable :: inum_interfaces_acoustic, inum_interfaces_elastic
+
+ double precision, dimension(:,:), allocatable :: buffer_send_faces_vector_acoustic
+ double precision, dimension(:,:), allocatable :: buffer_recv_faces_vector_acoustic
+ integer, dimension(:), allocatable :: tab_requests_send_recv_acoustic
+ double precision, dimension(:,:), allocatable :: buffer_send_faces_vector_elastic
+ double precision, dimension(:,:), allocatable :: buffer_recv_faces_vector_elastic
+ integer, dimension(:), allocatable :: tab_requests_send_recv_elastic
+ integer :: max_ibool_interfaces_size_acoustic, max_ibool_interfaces_size_elastic
+
+ integer, dimension(:,:), allocatable :: acoustic_surface
+ integer, dimension(:,:), allocatable :: acoustic_edges
+
+
+
+ integer :: ixmin, ixmax, izmin, izmax
+
+ integer :: ie, num_interface
+
+ integer :: nrecloc, irecloc
+ integer, dimension(:), allocatable :: recloc, which_proc_receiver
+
+ character(len=256) :: filename
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
!
!***********************************************************************
- open(IIN,file='OUTPUT_FILES/Database',status='old',action='read')
+#ifdef USE_MPI
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+#else
+ nproc = 1
+ myrank = 0
+
+#endif
+
+ write(prname,230)myrank
+230 format('./OUTPUT_FILES/Database',i5.5)
+ open(unit=IIN,file=prname,status='old',action='read')
+
+
! determine if we write to file instead of standard output
if(IOUT /= ISTANDARD_OUTPUT) open(IOUT,file='simulation_results.txt',status='unknown')
@@ -355,7 +425,7 @@
read(IIN,"(a80)") datlin
read(IIN,*) numat,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
read(IIN,"(a80)") datlin
- read(IIN,*) nelemabs,nelem_acoustic_surface
+ read(IIN,*) nelemabs,nelem_acoustic_surface,num_fluid_solid_edges
!
!---- allocate arrays
@@ -402,14 +472,14 @@
allocate(jend_right(nelemabs))
! --- allocate array for free surface condition in acoustic medium
- if(nelem_acoustic_surface <= 0) then
- nelem_acoustic_surface = 0
- allocate(ispecnum_acoustic_surface(1))
- allocate(iedgenum_acoustic_surface(1))
- else
- allocate(ispecnum_acoustic_surface(nelem_acoustic_surface))
- allocate(iedgenum_acoustic_surface(nelem_acoustic_surface))
- endif
+!!$ if(nelem_acoustic_surface <= 0) then
+!!$ nelem_acoustic_surface = 0
+!!$ allocate(ispecnum_acoustic_surface(1))
+!!$ allocate(iedgenum_acoustic_surface(1))
+!!$ else
+!!$ allocate(ispecnum_acoustic_surface(nelem_acoustic_surface))
+!!$ allocate(iedgenum_acoustic_surface(nelem_acoustic_surface))
+!!$ endif
!
!---- print element group main parameters
@@ -472,13 +542,52 @@
allocate(duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
+
!
+!---- read interfaces data
+!
+ print *, 'read the interfaces', myrank
+ read(IIN,"(a80)") datlin
+ read(IIN,*) ninterface, max_interface_size
+ if ( ninterface == 0 ) then
+ !allocate(my_neighbours(1))
+ !allocate(my_nelmnts_neighbours(1))
+ !allocate(my_interfaces(4,1,1))
+ !allocate(ibool_interfaces(NGLLX*1,1,1))
+ !allocate(nibool_interfaces(1,1))
+
+ else
+ allocate(my_neighbours(ninterface))
+ allocate(my_nelmnts_neighbours(ninterface))
+ allocate(my_interfaces(4,max_interface_size,ninterface))
+ allocate(ibool_interfaces_acoustic(NGLLX*max_interface_size,ninterface))
+ allocate(ibool_interfaces_elastic(NGLLX*max_interface_size,ninterface))
+ allocate(nibool_interfaces_acoustic(ninterface))
+ allocate(nibool_interfaces_elastic(ninterface))
+ allocate(inum_interfaces_acoustic(ninterface))
+ allocate(inum_interfaces_elastic(ninterface))
+
+ do num_interface = 1, ninterface
+ read(IIN,*) my_neighbours(num_interface), my_nelmnts_neighbours(num_interface)
+ do ie = 1, my_nelmnts_neighbours(num_interface)
+ read(IIN,*) my_interfaces(1,ie,num_interface), my_interfaces(2,ie,num_interface), &
+ my_interfaces(3,ie,num_interface), my_interfaces(4,ie,num_interface)
+
+ end do
+ end do
+ print *, 'end read the interfaces', myrank
+
+ end if
+
+
+!
!---- read absorbing boundary data
!
+ read(IIN,"(a80)") datlin
if(anyabs) then
- read(IIN,"(a80)") datlin
- do inum = 1,nelemabs
- read(IIN,*) numabsread,codeabsread(1),codeabsread(2),codeabsread(3),codeabsread(4)
+ do inum = 1,nelemabs
+ read(IIN,*) numabsread,codeabsread(1),codeabsread(2),codeabsread(3),codeabsread(4), ibegin_bottom(inum), iend_bottom(inum), &
+ jbegin_right(inum), jend_right(inum), ibegin_top(inum), iend_top(inum), jbegin_left(inum), jend_left(inum)
if(numabsread < 1 .or. numabsread > nspec) stop 'Wrong absorbing element number'
numabs(inum) = numabsread
codeabs(IBOTTOM,inum) = codeabsread(1)
@@ -493,20 +602,52 @@
!
!---- read acoustic free surface data
!
+ read(IIN,"(a80)") datlin
if(nelem_acoustic_surface > 0) then
- read(IIN,"(a80)") datlin
- do inum = 1,nelem_acoustic_surface
- read(IIN,*) numacoustread,iedgeacoustread
- if(numacoustread < 1 .or. numacoustread > nspec) stop 'Wrong acoustic free surface element number'
- if(iedgeacoustread < 1 .or. iedgeacoustread > NEDGES) stop 'Wrong acoustic free surface edge number'
- ispecnum_acoustic_surface(inum) = numacoustread
- iedgenum_acoustic_surface(inum) = iedgeacoustread
- enddo
+ allocate(acoustic_edges(4,nelem_acoustic_surface))
+ print *, 'POYU'
+ do inum = 1,nelem_acoustic_surface
+ read(IIN,*) acoustic_edges(1,inum), acoustic_edges(2,inum), acoustic_edges(3,inum), &
+ acoustic_edges(4,inum)
+ end do
+ print *, 'POYU'
+ allocate(acoustic_surface(5,nelem_acoustic_surface))
+ print *, 'POYU'
+ call construct_acoustic_surface ( nspec, ngnod, knods, nelem_acoustic_surface, &
+ acoustic_edges, acoustic_surface)
+ print *, 'POYU'
+!!$ read(IIN,*) numacoustread,iedgeacoustread
+!!$ if(numacoustread < 1 .or. numacoustread > nspec) stop 'Wrong acoustic free surface element number'
+!!$ if(iedgeacoustread < 1 .or. iedgeacoustread > NEDGES) stop 'Wrong acoustic free surface edge number'
+!!$ ispecnum_acoustic_surface(inum) = numacoustread
+!!$ iedgenum_acoustic_surface(inum) = iedgeacoustread
+!!$ enddo
write(IOUT,*)
write(IOUT,*) 'Number of free surface elements: ',nelem_acoustic_surface
endif
!
+!---- read acoustic elastic coupled edges
+!
+ read(IIN,"(a80)") datlin
+ if ( num_fluid_solid_edges > 0 ) then
+ allocate(fluid_solid_acoustic_ispec(num_fluid_solid_edges))
+ allocate(fluid_solid_acoustic_iedge(num_fluid_solid_edges))
+ allocate(fluid_solid_elastic_ispec(num_fluid_solid_edges))
+ allocate(fluid_solid_elastic_iedge(num_fluid_solid_edges))
+ do inum = 1, num_fluid_solid_edges
+ read(IIN,*) fluid_solid_acoustic_ispec(inum), fluid_solid_elastic_ispec(inum)
+ end do
+ else
+ allocate(fluid_solid_acoustic_ispec(1))
+ allocate(fluid_solid_acoustic_iedge(1))
+ allocate(fluid_solid_elastic_ispec(1))
+ allocate(fluid_solid_elastic_iedge(1))
+
+ end if
+
+
+!
!---- close input file
!
close(IIN)
@@ -560,13 +701,7 @@
if(nrec < 1) stop 'need at least one receiver'
-! allocate seismogram arrays
- allocate(sisux(NSTEP,nrec))
- allocate(sisuz(NSTEP,nrec))
-! to write seismograms in single precision SEP and double precision binary format
- allocate(buffer_binary_single(NSTEP*nrec))
- allocate(buffer_binary_double(NSTEP*nrec))
! receiver information
allocate(ispec_selected_rec(nrec))
@@ -576,6 +711,8 @@
allocate(gamma_receiver(nrec))
allocate(station_name(nrec))
allocate(network_name(nrec))
+ allocate(recloc(nrec))
+ allocate(which_proc_receiver(nrec))
! allocate 1-D Lagrange interpolators and derivatives
allocate(hxir(NGLLX))
@@ -706,25 +843,29 @@
if(source_type == 1) then
! collocated force source
call locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,source_type, &
- ix_source,iz_source,ispec_selected_source,iglob_source)
+ ix_source,iz_source,ispec_selected_source,iglob_source,is_proc_source,nb_proc_source)
! check that acoustic source is not exactly on the free surface because pressure is zero there
- do ispec_acoustic_surface = 1,nelem_acoustic_surface
- ispec = ispecnum_acoustic_surface(ispec_acoustic_surface)
- iedge = iedgenum_acoustic_surface(ispec_acoustic_surface)
- if(.not. elastic(ispec) .and. ispec == ispec_selected_source) then
- if((iedge == IBOTTOM .and. iz_source == 1) .or. &
- (iedge == ITOP .and. iz_source == NGLLZ) .or. &
- (iedge == ILEFT .and. ix_source == 1) .or. &
- (iedge == IRIGHT .and. ix_source == NGLLX)) &
- stop 'an acoustic source cannot be located exactly on the free surface because pressure is zero there'
- endif
- enddo
-
+ if ( is_proc_source == 1 ) then
+ do ispec_acoustic_surface = 1,nelem_acoustic_surface
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+ if( .not. elastic(ispec) .and. ispec == ispec_selected_source ) then
+ do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
+ do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
+ iglob = ibool(i,j,ispec)
+ if ( iglob_source == iglob ) then
+ stop 'an acoustic source cannot be located exactly on the free surface because pressure is zero there'
+ end if
+ end do
+ end do
+ endif
+ enddo
+ end if
+
else if(source_type == 2) then
! moment-tensor source
- call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
- ispec_selected_source,xi_source,gamma_source,coorg,knods,ngnod,npgeo)
+ call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
+ ispec_selected_source,is_proc_source,nb_proc_source,nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod,npgeo)
! compute source array for moment-tensor source
call compute_arrays_source(ispec_selected_source,xi_source,gamma_source,sourcearray, &
@@ -736,31 +877,52 @@
! locate receivers in the mesh
- call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,st_xval,st_zval,ispec_selected_rec, &
- xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo)
+ call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
+ st_xval,st_zval,ispec_selected_rec, &
+ xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo)
+! allocate seismogram arrays
+ allocate(sisux(NSTEP,nrecloc))
+ allocate(sisuz(NSTEP,nrecloc))
+
! check if acoustic receiver is exactly on the free surface because pressure is zero there
do ispec_acoustic_surface = 1,nelem_acoustic_surface
- ispec = ispecnum_acoustic_surface(ispec_acoustic_surface)
- iedge = iedgenum_acoustic_surface(ispec_acoustic_surface)
- do irec = 1,nrec
- if(.not. elastic(ispec) .and. ispec == ispec_selected_rec(irec)) then
- if((iedge == IBOTTOM .and. gamma_receiver(irec) < -0.99d0) .or. &
- (iedge == ITOP .and. gamma_receiver(irec) > 0.99d0) .or. &
- (iedge == ILEFT .and. xi_receiver(irec) < -0.99d0) .or. &
- (iedge == IRIGHT .and. xi_receiver(irec) > 0.99d0)) then
- if(seismotype == 4) then
- stop 'an acoustic pressure receiver cannot be located exactly on the free surface because pressure is zero there'
- else
- print *, '**********************************************************************'
- print *, '*** Warning: acoustic receiver located exactly on the free surface ***'
- print *, '*** Warning: tangential component will be zero there ***'
- print *, '**********************************************************************'
- print *
- endif
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+ ixmin = acoustic_surface(2,ispec_acoustic_surface)
+ ixmax = acoustic_surface(3,ispec_acoustic_surface)
+ izmin = acoustic_surface(4,ispec_acoustic_surface)
+ izmax = acoustic_surface(5,ispec_acoustic_surface)
+ do irecloc = 1,nrecloc
+ irec = recloc(irecloc)
+ if(.not. elastic(ispec) .and. ispec == ispec_selected_rec(irec)) then
+ if ( (izmin==1 .and. izmax==1 .and. ixmin==1 .and. ixmax==NGLLX .and. &
+ gamma_receiver(irec) < -0.99d0) .or.&
+ (izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==NGLLX .and. &
+ gamma_receiver(irec) > 0.99d0) .or.&
+ (izmin==1 .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==1 .and. &
+ xi_receiver(irec) < -0.99d0) .or.&
+ (izmin==1 .and. izmax==NGLLZ .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
+ xi_receiver(irec) > 0.99d0) .or.&
+ (izmin==1 .and. izmax==1 .and. ixmin==1 .and. ixmax==1 .and. &
+ gamma_receiver(irec) < -0.99d0 .and. xi_receiver(irec) < -0.99d0) .or.&
+ (izmin==1 .and. izmax==1 .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
+ gamma_receiver(irec) < -0.99d0 .and. xi_receiver(irec) > 0.99d0) .or.&
+ (izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==1 .and. &
+ gamma_receiver(irec) > 0.99d0 .and. xi_receiver(irec) < -0.99d0) .or.&
+ (izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
+ gamma_receiver(irec) > 0.99d0 .and. xi_receiver(irec) > 0.99d0) ) then
+ if(seismotype == 4) then
+ stop 'an acoustic pressure receiver cannot be located exactly on the free surface because pressure is zero there'
+ else
+ print *, '**********************************************************************'
+ print *, '*** Warning: acoustic receiver located exactly on the free surface ***'
+ print *, '*** Warning: tangential component will be zero there ***'
+ print *, '**********************************************************************'
+ print *
+ endif
+ endif
endif
- endif
- enddo
+ enddo
enddo
! define and store Lagrange interpolators at all the receivers
@@ -828,6 +990,59 @@
enddo
enddo
+#ifdef USE_MPI
+ if ( nproc > 1 ) then
+! preparing for MPI communications
+ call prepare_assemble_MPI (myrank,nspec,ibool, &
+ knods, ngnod, &
+ npoin, elastic, &
+ ninterface, max_interface_size, &
+ my_neighbours, my_nelmnts_neighbours, my_interfaces, &
+ ibool_interfaces_acoustic, ibool_interfaces_elastic, &
+ nibool_interfaces_acoustic, nibool_interfaces_elastic, &
+ inum_interfaces_acoustic, inum_interfaces_elastic, &
+ ninterface_acoustic, ninterface_elastic &
+ )
+
+ max_ibool_interfaces_size_acoustic = maxval(nibool_interfaces_acoustic(:))
+ max_ibool_interfaces_size_elastic = NDIM*maxval(nibool_interfaces_elastic(:))
+ allocate(tab_requests_send_recv_acoustic(ninterface_acoustic*2))
+ allocate(buffer_send_faces_vector_acoustic(max_ibool_interfaces_size_acoustic,ninterface_acoustic))
+ allocate(buffer_recv_faces_vector_acoustic(max_ibool_interfaces_size_acoustic,ninterface_acoustic))
+ allocate(tab_requests_send_recv_elastic(ninterface_elastic*2))
+ allocate(buffer_send_faces_vector_elastic(max_ibool_interfaces_size_elastic,ninterface_elastic))
+ allocate(buffer_recv_faces_vector_elastic(max_ibool_interfaces_size_elastic,ninterface_elastic))
+
+ call create_MPI_requests_SEND_RECV_acoustic(myrank, &
+ ninterface, ninterface_acoustic, &
+ nibool_interfaces_acoustic, &
+ my_neighbours, &
+ max_ibool_interfaces_size_acoustic, &
+ buffer_send_faces_vector_acoustic, &
+ buffer_recv_faces_vector_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ inum_interfaces_acoustic &
+ )
+
+ call create_MPI_requests_SEND_RECV_elastic(myrank, &
+ ninterface, ninterface_elastic, &
+ nibool_interfaces_elastic, &
+ my_neighbours, &
+ max_ibool_interfaces_size_elastic, &
+ buffer_send_faces_vector_elastic, &
+ buffer_recv_faces_vector_elastic, &
+ tab_requests_send_recv_elastic, &
+ inum_interfaces_elastic &
+ )
+
+ call assemble_MPI_scalar(myrank,rmass_inverse_acoustic, rmass_inverse_elastic,npoin, &
+ ninterface, max_interface_size, max_ibool_interfaces_size_acoustic, max_ibool_interfaces_size_elastic, &
+ ibool_interfaces_acoustic,ibool_interfaces_elastic, nibool_interfaces_acoustic,nibool_interfaces_elastic, my_neighbours)
+end if
+
+#endif
+
+
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
if(any_elastic) where(rmass_inverse_elastic <= 0.d0) rmass_inverse_elastic = 1.d0
if(any_acoustic) where(rmass_inverse_acoustic <= 0.d0) rmass_inverse_acoustic = 1.d0
@@ -839,11 +1054,13 @@
! check the mesh, stability and number of points per wavelength
call checkgrid(vpext,vsext,rhoext,density,elastcoef,ibool,kmato,coord,npoin,vpmin,vpmax, &
assign_external_model,nspec,numat,deltat,f0,t0,initialfield,time_function_type, &
- coorg,xinterp,zinterp,shape2D_display,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic)
+ coorg,xinterp,zinterp,shape2D_display,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic,myrank,nproc)
! convert receiver angle to radians
anglerec = anglerec * pi / 180.d0
+
+
!
!---- for color images
!
@@ -851,17 +1068,32 @@
if(output_color_image) then
! horizontal size of the image
- xmin_color_image = minval(coord(1,:))
- xmax_color_image = maxval(coord(1,:))
+ xmin_color_image_loc = minval(coord(1,:))
+ xmax_color_image_loc = maxval(coord(1,:))
! vertical size of the image, slightly increase it to go beyond maximum topography
- zmin_color_image = minval(coord(2,:))
- zmax_color_image = maxval(coord(2,:))
+ zmin_color_image_loc = minval(coord(2,:))
+ zmax_color_image_loc = maxval(coord(2,:))
+!
+ xmin_color_image = xmin_color_image_loc
+ xmax_color_image = xmax_color_image_loc
+ zmin_color_image = zmin_color_image_loc
+ zmax_color_image = zmax_color_image_loc
+
+#ifdef USE_MPI
+ call MPI_ALLREDUCE(xmin_color_image_loc, xmin_color_image, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE(xmax_color_image_loc, xmax_color_image, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE(zmin_color_image_loc, zmin_color_image, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE(zmax_color_image_loc, zmax_color_image, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+
+#endif
+
zmax_color_image = zmin_color_image + 1.05d0 * (zmax_color_image - zmin_color_image)
! compute number of pixels in the horizontal direction based on typical number
! of spectral elements in a given direction (may give bad results for very elongated models)
- NX_IMAGE_color = nint(sqrt(dble(npgeo))) * (NGLLX-1) + 1
+ !NX_IMAGE_color = nint(sqrt(dble(npgeo))) * (NGLLX-1) + 1
+ NX_IMAGE_color = 800
! compute number of pixels in the vertical direction based on ratio of sizes
NZ_IMAGE_color = nint(NX_IMAGE_color * (zmax_color_image - zmin_color_image) / (xmax_color_image - xmin_color_image))
@@ -886,92 +1118,191 @@
iglob_image_color(:,:) = -1
-! loop on all the grid points to map them to a pixel in the image
- do n=1,npoin
+ !!!!!!
+ nb_pixel_loc = 0
+ do ispec = 1, nspec
+
+ do k = 1, 4
+ elmnt_coords(1,k) = coorg(1,knods(k,ispec))
+ elmnt_coords(2,k) = coorg(2,knods(k,ispec))
-! compute the coordinates of this pixel
- i = nint((coord(1,n) - xmin_color_image) / size_pixel_horizontal + 1)
- j = nint((coord(2,n) - zmin_color_image) / size_pixel_vertical + 1)
+ end do
+
+ min_i = floor(minval((elmnt_coords(1,:) - xmin_color_image))/size_pixel_horizontal) + 1
+ max_i = ceiling(maxval((elmnt_coords(1,:) - xmin_color_image))/size_pixel_horizontal) + 1
+ min_j = floor(minval((elmnt_coords(2,:) - zmin_color_image))/size_pixel_vertical) + 1
+ max_j = ceiling(maxval((elmnt_coords(2,:) - zmin_color_image))/size_pixel_vertical) + 1
+
+ do j = min_j, max_j
+ do i = min_i, max_i
+ i_coord = (i-1)*size_pixel_horizontal + xmin_color_image
+ j_coord = (j-1)*size_pixel_vertical + zmin_color_image
+
+ call is_in_convex_quadrilateral( elmnt_coords, i_coord, j_coord, pixel_is_in)
-! avoid edge effects
- if(i < 1) i = 1
- if(i > NX_IMAGE_color) i = NX_IMAGE_color
+ if ( pixel_is_in ) then
+ dist_min_pixel = HUGEVAL
+ do k = 1, NGLLX
+ do l = 1, NGLLZ
+ iglob = ibool(k,l,ispec)
+ dist_pixel = (coord(1,iglob)-i_coord)**2 + (coord(2,iglob)-j_coord)**2
+ if (dist_pixel < dist_min_pixel) then
+ dist_min_pixel = dist_pixel
+ iglob_image_color(i,j) = iglob
+
+ end if
+
+ end do
+ end do
+ if ( dist_min_pixel >= HUGEVAL ) then
+ stop 'Error in detecting pixel for color image'
+
+ end if
+ nb_pixel_loc = nb_pixel_loc + 1
+
+ end if
+
+ end do
+ end do
+ end do
- if(j < 1) j = 1
- if(j > NZ_IMAGE_color) j = NZ_IMAGE_color
+ print *, 'POYOP', nb_pixel_loc, myrank
+ allocate(num_pixel_loc(nb_pixel_loc))
+
+ nb_pixel_loc = 0
+ do i = 1, NX_IMAGE_color
+ do j = 1, NZ_IMAGE_color
+ if ( iglob_image_color(i,j) /= -1 ) then
+ nb_pixel_loc = nb_pixel_loc + 1
+ num_pixel_loc(nb_pixel_loc) = (j-1)*NX_IMAGE_color + i
+
+ end if
+
+ end do
+ end do
+
+#ifdef USE_MPI
+ allocate(nb_pixel_per_proc(nproc))
+
+ call MPI_GATHER( nb_pixel_loc, 1, MPI_INTEGER, nb_pixel_per_proc(1), 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ier)
+
+ if ( myrank == 0 ) then
+ allocate(num_pixel_recv(maxval(nb_pixel_per_proc(:)),nproc))
+ allocate(data_pixel_recv(maxval(nb_pixel_per_proc(:))))
+ end if
-! assign this point to this pixel
- iglob_image_color(i,j) = n
+ allocate(data_pixel_send(nb_pixel_loc))
+ if ( nproc > 1 ) then
+ if ( myrank == 0 ) then
+
+ do iproc = 1, nproc-1
- enddo
+ call MPI_RECV(num_pixel_recv(1,iproc+1),nb_pixel_per_proc(iproc+1), MPI_INTEGER, &
+ iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
+ do k = 1, nb_pixel_per_proc(iproc+1)
+ j = ceiling(real(num_pixel_recv(k,iproc+1)) / real(NX_IMAGE_color))
+ i = num_pixel_recv(k,iproc+1) - (j-1)*NX_IMAGE_color
+ iglob_image_color(i,j) = iproc
+
+ end do
+ end do
+
+ else
+ call MPI_SEND(num_pixel_loc(1),nb_pixel_loc,MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
+
+ end if
-! locate missing pixels based on a minimum distance criterion
-! cannot do more than two iterations typically because some pixels must never be found
-! because they do not exist (for instance if they are located above topography)
- do count_passes = 1,2
+ end if
+#endif
+
+
- print *,'pass ',count_passes,' to locate the missing pixels of color images'
+!!$! loop on all the grid points to map them to a pixel in the image
+!!$ do n=1,npoin
+!!$
+!!$! compute the coordinates of this pixel
+!!$ i = nint((coord(1,n) - xmin_color_image) / size_pixel_horizontal + 1)
+!!$ j = nint((coord(2,n) - zmin_color_image) / size_pixel_vertical + 1)
+!!$
+!!$! avoid edge effects
+!!$ if(i < 1) i = 1
+!!$ if(i > NX_IMAGE_color) i = NX_IMAGE_color
+!!$
+!!$ if(j < 1) j = 1
+!!$ if(j > NZ_IMAGE_color) j = NZ_IMAGE_color
+!!$
+!!$! assign this point to this pixel
+!!$ iglob_image_color(i,j) = n
+!!$
+!!$ enddo
+!!$
+!!$! locate missing pixels based on a minimum distance criterion
+!!$! cannot do more than two iterations typically because some pixels must never be found
+!!$! because they do not exist (for instance if they are located above topography)
+!!$ do count_passes = 1,20
+!!$
+!!$ print *,'pass ',count_passes,' to locate the missing pixels of color images'
+!!$
+!!$ copy_iglob_image_color(:,:) = iglob_image_color(:,:)
+!!$
+!!$ do j = 1,NZ_IMAGE_color
+!!$ do i = 1,NX_IMAGE_color
+!!$
+!!$ if(copy_iglob_image_color(i,j) == -1) then
+!!$
+!!$ iplus1 = i + 1
+!!$ iminus1 = i - 1
+!!$
+!!$ jplus1 = j + 1
+!!$ jminus1 = j - 1
+!!$
+!!$! avoid edge effects
+!!$ if(iminus1 < 1) iminus1 = 1
+!!$ if(iplus1 > NX_IMAGE_color) iplus1 = NX_IMAGE_color
+!!$
+!!$ if(jminus1 < 1) jminus1 = 1
+!!$ if(jplus1 > NZ_IMAGE_color) jplus1 = NZ_IMAGE_color
+!!$
+!!$! use neighbors of this pixel to fill the holes
+!!$
+!!$! horizontal
+!!$ if(copy_iglob_image_color(iplus1,j) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(iplus1,j)
+!!$
+!!$ else if(copy_iglob_image_color(iminus1,j) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(iminus1,j)
+!!$
+!!$! vertical
+!!$ else if(copy_iglob_image_color(i,jplus1) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(i,jplus1)
+!!$
+!!$ else if(copy_iglob_image_color(i,jminus1) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(i,jminus1)
+!!$
+!!$! diagonal
+!!$ else if(copy_iglob_image_color(iminus1,jminus1) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(iminus1,jminus1)
+!!$
+!!$ else if(copy_iglob_image_color(iplus1,jminus1) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(iplus1,jminus1)
+!!$
+!!$ else if(copy_iglob_image_color(iminus1,jplus1) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(iminus1,jplus1)
+!!$
+!!$ else if(copy_iglob_image_color(iplus1,jplus1) /= -1) then
+!!$ iglob_image_color(i,j) = copy_iglob_image_color(iplus1,jplus1)
+!!$
+!!$ endif
+!!$
+!!$ endif
+!!$
+!!$ enddo
+!!$ enddo
+!!$
+!!$ enddo
+!!$
+!!$ deallocate(copy_iglob_image_color)
- copy_iglob_image_color(:,:) = iglob_image_color(:,:)
-
- do j = 1,NZ_IMAGE_color
- do i = 1,NX_IMAGE_color
-
- if(copy_iglob_image_color(i,j) == -1) then
-
- iplus1 = i + 1
- iminus1 = i - 1
-
- jplus1 = j + 1
- jminus1 = j - 1
-
-! avoid edge effects
- if(iminus1 < 1) iminus1 = 1
- if(iplus1 > NX_IMAGE_color) iplus1 = NX_IMAGE_color
-
- if(jminus1 < 1) jminus1 = 1
- if(jplus1 > NZ_IMAGE_color) jplus1 = NZ_IMAGE_color
-
-! use neighbors of this pixel to fill the holes
-
-! horizontal
- if(copy_iglob_image_color(iplus1,j) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(iplus1,j)
-
- else if(copy_iglob_image_color(iminus1,j) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(iminus1,j)
-
-! vertical
- else if(copy_iglob_image_color(i,jplus1) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(i,jplus1)
-
- else if(copy_iglob_image_color(i,jminus1) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(i,jminus1)
-
-! diagonal
- else if(copy_iglob_image_color(iminus1,jminus1) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(iminus1,jminus1)
-
- else if(copy_iglob_image_color(iplus1,jminus1) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(iplus1,jminus1)
-
- else if(copy_iglob_image_color(iminus1,jplus1) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(iminus1,jplus1)
-
- else if(copy_iglob_image_color(iplus1,jplus1) /= -1) then
- iglob_image_color(i,j) = copy_iglob_image_color(iplus1,jplus1)
-
- endif
-
- endif
-
- enddo
- enddo
-
- enddo
-
- deallocate(copy_iglob_image_color)
-
write(IOUT,*) 'done locating all the pixels of color images'
endif
@@ -1071,31 +1402,33 @@
! output absolute time in third column, in case user wants to check it as well
write(55,*) sngl(time),sngl(source_time_function(it)),sngl(time-t0)
- enddo
+ enddo
close(55)
+
+ source_time_function(:) = source_time_function(:) / nb_proc_source
else
allocate(source_time_function(1))
- endif
+ endif
!
!---- check that no element has both acoustic free surface and top absorbing surface
!
- do ispec_acoustic_surface = 1,nelem_acoustic_surface
- ispec = ispecnum_acoustic_surface(ispec_acoustic_surface)
- iedge = iedgenum_acoustic_surface(ispec_acoustic_surface)
- if(elastic(ispec)) then
- stop 'elastic element detected in acoustic free surface'
- else
- do inum = 1,nelemabs
- if(numabs(inum) == ispec .and. codeabs(iedge,inum)) &
- stop 'acoustic free surface cannot be both absorbing and free'
- enddo
- endif
- enddo
+!!$ do ispec_acoustic_surface = 1,nelem_acoustic_surface
+!!$ ispec = ispecnum_acoustic_surface(ispec_acoustic_surface)
+!!$ iedge = iedgenum_acoustic_surface(ispec_acoustic_surface)
+!!$ if(elastic(ispec)) then
+!!$ stop 'elastic element detected in acoustic free surface'
+!!$ else
+!!$ do inum = 1,nelemabs
+!!$ if(numabs(inum) == ispec .and. codeabs(iedge,inum)) &
+!!$ stop 'acoustic free surface cannot be both absorbing and free'
+!!$ enddo
+!!$ endif
+!!$ enddo
! determine if coupled fluid-solid simulation
coupled_acoustic_elastic = any_acoustic .and. any_elastic
@@ -1104,8 +1437,6 @@
! very basic algorithm in O(nspec^2), simple double loop on the elements
! and then loop on the four corners of each of the two elements, could be signficantly improved
- num_fluid_solid_edges_alloc = 0
-
if(coupled_acoustic_elastic) then
print *
print *,'Mixed acoustic/elastic simulation'
@@ -1158,56 +1489,11 @@
enddo
-! double loop on all the elements
- do ispec_acoustic = 1, nspec
- do ispec_elastic = 1, nspec
-! one element must be acoustic and the other must be elastic
-! use acoustic element as master to avoid double detection of the same pair (one on each side)
- if(ispec_acoustic /= ispec_elastic .and. .not. elastic(ispec_acoustic) .and. elastic(ispec_elastic)) then
+ do inum = 1, num_fluid_solid_edges
+ ispec_acoustic = fluid_solid_acoustic_ispec(inum)
+ ispec_elastic = fluid_solid_elastic_ispec(inum)
-! loop on the four edges of the two elements
- do iedge_acoustic = 1,NEDGES
- do iedge_elastic = 1,NEDGES
-
-! error if the two edges match in direct order
- if(ibool(i_begin(iedge_acoustic),j_begin(iedge_acoustic),ispec_acoustic) == &
- ibool(i_begin(iedge_elastic),j_begin(iedge_elastic),ispec_elastic) .and. &
- ibool(i_end(iedge_acoustic),j_end(iedge_acoustic),ispec_acoustic) == &
- ibool(i_end(iedge_elastic),j_end(iedge_elastic),ispec_elastic)) &
- stop 'topology error (non-inverted coupled elements) found in fluid/solid edge detection'
-
-! the two edges can match in inverse order
- if(ibool(i_begin(iedge_acoustic),j_begin(iedge_acoustic),ispec_acoustic) == &
- ibool(i_end(iedge_elastic),j_end(iedge_elastic),ispec_elastic) .and. &
- ibool(i_end(iedge_acoustic),j_end(iedge_acoustic),ispec_acoustic) == &
- ibool(i_begin(iedge_elastic),j_begin(iedge_elastic),ispec_elastic)) &
- num_fluid_solid_edges_alloc = num_fluid_solid_edges_alloc + 1
-
- enddo
- enddo
-
- endif
-
- enddo
- enddo
-
- print *,'Number of fluid/solid edges detected in mesh = ',num_fluid_solid_edges_alloc
-
-! allocate arrays for fluid/solid matching
- allocate(fluid_solid_acoustic_ispec(num_fluid_solid_edges_alloc))
- allocate(fluid_solid_acoustic_iedge(num_fluid_solid_edges_alloc))
- allocate(fluid_solid_elastic_ispec(num_fluid_solid_edges_alloc))
- allocate(fluid_solid_elastic_iedge(num_fluid_solid_edges_alloc))
-
-! double loop on all the elements
- print *,'Creating fluid/solid edge topology...'
-
- num_fluid_solid_edges = 0
-
- do ispec_acoustic = 1, nspec
- do ispec_elastic = 1, nspec
-
! one element must be acoustic and the other must be elastic
! use acoustic element as master to avoid double detection of the same pair (one on each side)
if(ispec_acoustic /= ispec_elastic .and. .not. elastic(ispec_acoustic) .and. elastic(ispec_elastic)) then
@@ -1218,27 +1504,24 @@
! store the matching topology if the two edges match in inverse order
if(ibool(i_begin(iedge_acoustic),j_begin(iedge_acoustic),ispec_acoustic) == &
- ibool(i_end(iedge_elastic),j_end(iedge_elastic),ispec_elastic) .and. &
- ibool(i_end(iedge_acoustic),j_end(iedge_acoustic),ispec_acoustic) == &
- ibool(i_begin(iedge_elastic),j_begin(iedge_elastic),ispec_elastic)) then
- num_fluid_solid_edges = num_fluid_solid_edges + 1
- fluid_solid_acoustic_ispec(num_fluid_solid_edges) = ispec_acoustic
- fluid_solid_acoustic_iedge(num_fluid_solid_edges) = iedge_acoustic
- fluid_solid_elastic_ispec(num_fluid_solid_edges) = ispec_elastic
- fluid_solid_elastic_iedge(num_fluid_solid_edges) = iedge_elastic
+ ibool(i_end(iedge_elastic),j_end(iedge_elastic),ispec_elastic) .and. &
+ ibool(i_end(iedge_acoustic),j_end(iedge_acoustic),ispec_acoustic) == &
+ ibool(i_begin(iedge_elastic),j_begin(iedge_elastic),ispec_elastic)) then
+ fluid_solid_acoustic_iedge(inum) = iedge_acoustic
+ fluid_solid_elastic_iedge(inum) = iedge_elastic
! print *,'edge ',iedge_acoustic,' of acoustic element ',ispec_acoustic, &
! ' is in contact with edge ',iedge_elastic,' of elastic element ',ispec_elastic
- endif
+ endif
- enddo
+ enddo
enddo
- endif
+ endif
- enddo
enddo
+
- if(num_fluid_solid_edges /= num_fluid_solid_edges_alloc) stop 'error in creation of arrays for fluid/solid matching'
+ !if(num_fluid_solid_edges /= num_fluid_solid_edges_alloc) stop 'error in creation of arrays for fluid/solid matching'
! make sure fluid/solid matching has been perfectly detected: check that the grid points
! have the same physical coordinates
@@ -1282,27 +1565,23 @@
else
-! allocate dummy arrays for fluid/solid matching if purely acoustic or purely elastic
- allocate(fluid_solid_acoustic_ispec(1))
- allocate(fluid_solid_acoustic_iedge(1))
- allocate(fluid_solid_elastic_ispec(1))
- allocate(fluid_solid_elastic_iedge(1))
+
endif
-
+
! default values for acoustic absorbing edges
- ibegin_bottom(:) = 1
- ibegin_top(:) = 1
+!!$ ibegin_bottom(:) = 1
+!!$ ibegin_top(:) = 1
+!!$
+!!$ iend_bottom(:) = NGLLX
+!!$ iend_top(:) = NGLLX
+!!$
+!!$ jbegin_left(:) = 1
+!!$ jbegin_right(:) = 1
+!!$
+!!$ jend_left(:) = NGLLZ
+!!$ jend_right(:) = NGLLZ
- iend_bottom(:) = NGLLX
- iend_top(:) = NGLLX
-
- jbegin_left(:) = 1
- jbegin_right(:) = 1
-
- jend_left(:) = NGLLZ
- jend_right(:) = NGLLZ
-
! exclude common points between acoustic absorbing edges and acoustic/elastic matching interfaces
if(coupled_acoustic_elastic .and. anyabs) then
@@ -1322,8 +1601,8 @@
! if acoustic absorbing element and acoustic/elastic coupled element is the same
if(ispec_acoustic == ispec) then
-
- if(iedge_acoustic == IBOTTOM) then
+
+ if(iedge_acoustic == IBOTTOM) then
jbegin_left(ispecabs) = 2
jbegin_right(ispecabs) = 2
endif
@@ -1377,7 +1656,7 @@
! *********************************************************
do it = 1,NSTEP
-
+
! compute current time
time = (it-1)*deltat
@@ -1395,16 +1674,16 @@
potential_dot_dot_acoustic = ZERO
! free surface for an acoustic medium
- call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,ispecnum_acoustic_surface,iedgenum_acoustic_surface, &
- ibool,nelem_acoustic_surface,npoin,nspec)
+ call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
+ potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
! *********************************************************
! ************* compute forces for the acoustic elements
! *********************************************************
call compute_forces_acoustic(npoin,nspec,nelemabs,numat, &
- iglob_source,ispec_selected_source,source_type,it,NSTEP,anyabs, &
+ iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs, &
assign_external_model,initialfield,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
@@ -1413,13 +1692,13 @@
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right)
- endif ! end of test if any acoustic element
+ endif ! end of test if any acoustic element
! *********************************************************
! ************* add coupling with the elastic side
! *********************************************************
- if(coupled_acoustic_elastic) then
+ if(coupled_acoustic_elastic) then
! loop on all the coupling edges
do inum = 1,num_fluid_solid_edges
@@ -1479,8 +1758,40 @@
enddo
- endif
+ endif
+!!$ potential_dot_dot_acoustic = 1
+!!$ write(filename,282)myrank
+!!$282 format('OUTPUT_FILES/acoustic_before',i5.5)
+!!$ call gnuplot_array ( npoin, nspec, filename, potential_dot_dot_acoustic, coord, ibool)
+
+#ifdef USE_MPI
+ if ( nproc > 1 .and. any_acoustic .and. ninterface_acoustic > 0) then
+ call assemble_MPI_vector_acoustic_start(myrank,potential_dot_dot_acoustic,npoin, &
+ ninterface, ninterface_acoustic, &
+ inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_acoustic,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ buffer_send_faces_vector_acoustic, &
+ buffer_recv_faces_vector_acoustic &
+ )
+ call assemble_MPI_vector_acoustic_wait(myrank,potential_dot_dot_acoustic,npoin, &
+ ninterface, ninterface_acoustic, &
+ inum_interfaces_acoustic, &
+ max_interface_size, max_ibool_interfaces_size_acoustic,&
+ ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
+ tab_requests_send_recv_acoustic, &
+ buffer_send_faces_vector_acoustic, &
+ buffer_recv_faces_vector_acoustic &
+ )
+ end if
+#endif
+!!$ write(filename,283)myrank
+!!$283 format('OUTPUT_FILES/acoustic_after',i5.5)
+!!$ call gnuplot_array ( npoin, nspec, filename, potential_dot_dot_acoustic, coord, ibool)
+
+
! ************************************************************************************
! ************* multiply by the inverse of the mass matrix and update velocity
! ************************************************************************************
@@ -1492,17 +1803,17 @@
! free surface for an acoustic medium
call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,ispecnum_acoustic_surface,iedgenum_acoustic_surface, &
+ potential_acoustic,acoustic_surface, &
ibool,nelem_acoustic_surface,npoin,nspec)
endif
-
+
! *********************************************************
! ************* main solver for the elastic elements
! *********************************************************
- if(any_elastic) &
+ if(any_elastic) &
call compute_forces_elastic(npoin,nspec,nelemabs,numat,iglob_source, &
- ispec_selected_source,source_type,it,NSTEP,anyabs,assign_external_model, &
+ ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic,density,elastcoef,xix,xiz,gammax,gammaz, &
@@ -1582,16 +1893,55 @@
endif
+
+!!$ accel_elastic = 1
+!!$ write(filename,280)myrank
+!!$280 format('OUTPUT_FILES/accel_elastic_before',i5.5)
+!!$ call gnuplot_array ( npoin, nspec, filename, accel_elastic(2,:), coord, ibool)
+
! ************************************************************************************
! ************* multiply by the inverse of the mass matrix and update velocity
! ************************************************************************************
+#ifdef USE_MPI
+ if ( nproc > 1 .and. any_elastic .and. ninterface_elastic > 0) then
+ call assemble_MPI_vector_elastic_start(myrank,accel_elastic,npoin, &
+ ninterface, ninterface_elastic, &
+ inum_interfaces_elastic, &
+ max_interface_size, max_ibool_interfaces_size_elastic,&
+ ibool_interfaces_elastic, nibool_interfaces_elastic, &
+ tab_requests_send_recv_elastic, &
+ buffer_send_faces_vector_elastic, &
+ buffer_recv_faces_vector_elastic &
+ )
+ call assemble_MPI_vector_elastic_wait(myrank,accel_elastic,npoin, &
+ ninterface, ninterface_elastic, &
+ inum_interfaces_elastic, &
+ max_interface_size, max_ibool_interfaces_size_elastic,&
+ ibool_interfaces_elastic, nibool_interfaces_elastic, &
+ tab_requests_send_recv_elastic, &
+ buffer_send_faces_vector_elastic, &
+ buffer_recv_faces_vector_elastic &
+ )
+ end if
+#endif
+
+!!$ write(filename,281)myrank
+!!$281 format('OUTPUT_FILES/accel_elastic_after',i5.5)
+!!$ call gnuplot_array ( npoin, nspec, filename, accel_elastic(2,:), coord, ibool)
+!!$#ifdef USE_MPI
+!!$ call MPI_BARRIER(MPI_COMM_WORLD, ier)
+!!$#endif
+!!$
+!!$ stop
+
if(any_elastic) then
accel_elastic(1,:) = accel_elastic(1,:) * rmass_inverse_elastic
accel_elastic(2,:) = accel_elastic(2,:) * rmass_inverse_elastic
veloc_elastic = veloc_elastic + deltatover2*accel_elastic
endif
+
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
@@ -1620,8 +1970,10 @@
endif
! loop on all the receivers to compute and store the seismograms
- do irec = 1,nrec
-
+ do irecloc = 1,nrecloc
+
+ irec = recloc(irecloc)
+
ispec = ispec_selected_rec(irec)
! compute pressure in this element if needed
@@ -1695,11 +2047,11 @@
! rotate seismogram components if needed, except if recording pressure, which is a scalar
if(seismotype /= 4) then
- sisux(it,irec) = cosrot*valux + sinrot*valuz
- sisuz(it,irec) = - sinrot*valux + cosrot*valuz
+ sisux(it,irecloc) = cosrot*valux + sinrot*valuz
+ sisuz(it,irecloc) = - sinrot*valux + cosrot*valuz
else
- sisux(it,irec) = valux
- sisuz(it,irec) = ZERO
+ sisux(it,irecloc) = valux
+ sisuz(it,irecloc) = ZERO
endif
enddo
@@ -1723,14 +2075,14 @@
call compute_vector_whole_medium(potential_acoustic,displ_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
- call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
- it,deltat,coorg,xinterp,zinterp,shape2D_display, &
- Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
- numabs,codeabs,anyabs,simulation_title,npoin,npgeo,vpmin,vpmax,nrec, &
- colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
- boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
- nspec,ngnod,coupled_acoustic_elastic,any_acoustic,plot_lowerleft_corner_only, &
- fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges)
+!!$ call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
+!!$ it,deltat,coorg,xinterp,zinterp,shape2D_display, &
+!!$ Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
+!!$ numabs,codeabs,anyabs,simulation_title,npoin,npgeo,vpmin,vpmax,nrec, &
+!!$ colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
+!!$ boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
+!!$ nspec,ngnod,coupled_acoustic_elastic,any_acoustic,plot_lowerleft_corner_only, &
+!!$ fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges)
else if(imagetype == 2) then
@@ -1739,14 +2091,14 @@
call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
- call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
- it,deltat,coorg,xinterp,zinterp,shape2D_display, &
- Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
- numabs,codeabs,anyabs,simulation_title,npoin,npgeo,vpmin,vpmax,nrec, &
- colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
- boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
- nspec,ngnod,coupled_acoustic_elastic,any_acoustic,plot_lowerleft_corner_only, &
- fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges)
+!!$ call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
+!!$ it,deltat,coorg,xinterp,zinterp,shape2D_display, &
+!!$ Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
+!!$ numabs,codeabs,anyabs,simulation_title,npoin,npgeo,vpmin,vpmax,nrec, &
+!!$ colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
+!!$ boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
+!!$ nspec,ngnod,coupled_acoustic_elastic,any_acoustic,plot_lowerleft_corner_only, &
+!!$ fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges)
else if(imagetype == 3) then
@@ -1755,14 +2107,14 @@
call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
- call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
- it,deltat,coorg,xinterp,zinterp,shape2D_display, &
- Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
- numabs,codeabs,anyabs,simulation_title,npoin,npgeo,vpmin,vpmax,nrec, &
- colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
- boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
- nspec,ngnod,coupled_acoustic_elastic,any_acoustic,plot_lowerleft_corner_only, &
- fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges)
+!!$ call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
+!!$ it,deltat,coorg,xinterp,zinterp,shape2D_display, &
+!!$ Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
+!!$ numabs,codeabs,anyabs,simulation_title,npoin,npgeo,vpmin,vpmax,nrec, &
+!!$ colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
+!!$ boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
+!!$ nspec,ngnod,coupled_acoustic_elastic,any_acoustic,plot_lowerleft_corner_only, &
+!!$ fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges)
else if(imagetype == 4) then
@@ -1818,25 +2170,71 @@
endif
image_color_data(:,:) = 0.d0
- do j = 1,NZ_IMAGE_color
- do i = 1,NX_IMAGE_color
- iglob = iglob_image_color(i,j)
-! draw vertical component of field
-! or pressure which is stored in the same array used as temporary storage
- if(iglob /= -1) image_color_data(i,j) = vector_field_display(2,iglob)
- enddo
- enddo
+
+ do k = 1, nb_pixel_loc
+ j = ceiling(real(num_pixel_loc(k)) / real(NX_IMAGE_color))
+ i = num_pixel_loc(k) - (j-1)*NX_IMAGE_color
+ image_color_data(i,j) = vector_field_display(2,iglob_image_color(i,j))
+
+ end do
+
+#ifdef USE_MPI
+ if ( nproc > 1 ) then
+ if ( myrank == 0 ) then
+
+ do iproc = 1, nproc-1
- call create_color_image(image_color_data,iglob_image_color,NX_IMAGE_color,NZ_IMAGE_color,it,cutsnaps)
+ call MPI_RECV(data_pixel_recv(1),nb_pixel_per_proc(iproc+1), MPI_DOUBLE_PRECISION, &
+ iproc, 43, MPI_COMM_WORLD, request_mpi_status, ier)
+
+ do k = 1, nb_pixel_per_proc(iproc+1)
+ j = ceiling(real(num_pixel_recv(k,iproc+1)) / real(NX_IMAGE_color))
+ i = num_pixel_recv(k,iproc+1) - (j-1)*NX_IMAGE_color
+ image_color_data(i,j) = data_pixel_recv(k)
+
+ end do
+ end do
+
+ else
+ do k = 1, nb_pixel_loc
+ j = ceiling(real(num_pixel_loc(k)) / real(NX_IMAGE_color))
+ i = num_pixel_loc(k) - (j-1)*NX_IMAGE_color
+ data_pixel_send(k) = vector_field_display(2,iglob_image_color(i,j))
+
+ end do
+
+ call MPI_SEND(data_pixel_send(1),nb_pixel_loc,MPI_DOUBLE_PRECISION, 0, 43, MPI_COMM_WORLD, ier)
+
+ end if
+ end if
+! call MPI_BARRIER(MPI_COMM_WORLD, ier)
- write(IOUT,*) 'Color image created'
+#endif
+
+!!$ do j = 1,NZ_IMAGE_color
+!!$ do i = 1,NX_IMAGE_color
+!!$ iglob = iglob_image_color(i,j)
+!!$! draw vertical component of field
+!!$! or pressure which is stored in the same array used as temporary storage
+!!$ if(iglob /= -1) image_color_data(i,j) = vector_field_display(2,iglob)
+!!$ enddo
+!!$ enddo
+ if ( myrank == 0 ) then
+ call create_color_image(image_color_data,iglob_image_color,NX_IMAGE_color,NZ_IMAGE_color,it,cutsnaps)
+
+ write(IOUT,*) 'Color image created'
+
+ end if
+
endif
-
+
!---- save temporary or final seismograms
- call write_seismograms(sisux,sisuz,station_name,network_name,NSTEP, &
- nrec,deltat,seismotype,st_xval,it,t0,buffer_binary_single,buffer_binary_double)
-
+ if ( it == NSTEP ) then
+ call write_seismograms(sisux,sisuz,station_name,network_name,NSTEP, &
+ nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0)
+ end if
+
! count elapsed wall-clock time
call date_and_time(datein,timein,zone,time_values)
! time_values(3): day of the month
@@ -1942,3 +2340,68 @@
end program specfem2D
+
+
+
+subroutine gnuplot_array ( npoin, nspec, filename, array, coord, ibool)
+
+ implicit none
+ include "constants.h"
+
+ integer :: npoin
+ integer :: nspec
+ character(len=256) :: filename
+ double precision, dimension(npoin) :: array
+ double precision, dimension(NDIM,npoin) :: coord
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+
+ integer :: i, j, ispec
+ integer :: ipoint
+
+ open(unit=90, file=trim(filename), status='unknown', action='write')
+ do ipoint = 1, npoin
+
+ write(90,*) coord(1,ipoint), coord(2,ipoint), array(ipoint)
+ end do
+ close(90)
+
+
+end subroutine gnuplot_array
+
+
+subroutine is_in_convex_quadrilateral ( elmnt_coords, x_coord, z_coord, is_in)
+
+ implicit none
+
+ double precision, dimension(2,4) :: elmnt_coords
+ double precision, intent(in) :: x_coord, z_coord
+ logical, intent(out) :: is_in
+
+ real :: x1, x2, x3, x4, z1, z2, z3, z4
+ real :: normal1, normal2, normal3, normal4
+
+
+ x1 = elmnt_coords(1,1)
+ x2 = elmnt_coords(1,2)
+ x3 = elmnt_coords(1,3)
+ x4 = elmnt_coords(1,4)
+ z1 = elmnt_coords(2,1)
+ z2 = elmnt_coords(2,2)
+ z3 = elmnt_coords(2,3)
+ z4 = elmnt_coords(2,4)
+
+ normal1 = (z_coord-z1) * (x2-x1) - (x_coord-x1) * (z2-z1)
+ normal2 = (z_coord-z2) * (x3-x2) - (x_coord-x2) * (z3-z2)
+ normal3 = (z_coord-z3) * (x4-x3) - (x_coord-x3) * (z4-z3)
+ normal4 = (z_coord-z4) * (x1-x4) - (x_coord-x4) * (z1-z4)
+
+ if ( (normal1 < 0) .or. (normal2 < 0) .or. (normal3 < 0) .or. (normal4 < 0) ) then
+ is_in = .false.
+ !print *, 'normal', normal1, normal2, normal3, normal4
+ else
+ is_in = .true.
+ end if
+
+
+
+end subroutine is_in_convex_quadrilateral
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-06-26 00:36:41 UTC (rev 8523)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-12-07 23:53:26 UTC (rev 8524)
@@ -14,16 +14,22 @@
! write seismograms to text files
subroutine write_seismograms(sisux,sisuz,station_name,network_name, &
- NSTEP,nrec,deltat,seismotype,st_xval,it,t0,buffer_binary_single,buffer_binary_double)
+ NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0)
implicit none
include "constants.h"
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
integer nrec,NSTEP,it,seismotype
double precision t0,deltat
+
+ integer, intent(in) :: nrecloc,myrank
+ integer, dimension(nrec),intent(in) :: which_proc_receiver
- double precision, dimension(NSTEP,nrec) :: sisux,sisuz
+ double precision, dimension(NSTEP,nrecloc), intent(in) :: sisux,sisuz
double precision st_xval(nrec)
@@ -37,12 +43,18 @@
character(len=150) sisname
! to write seismograms in single precision SEP and double precision binary format
- real(kind=4), dimension(NSTEP*nrec) :: buffer_binary_single
- double precision, dimension(NSTEP*nrec) :: buffer_binary_double
+ double precision, dimension(:,:), allocatable :: buffer_binary
! scaling factor for Seismic Unix xsu dislay
double precision, parameter :: FACTORXSU = 1.d0
+
+ integer :: irecloc
+ integer :: ierror
+#ifdef USE_MPI
+ integer, dimension(MPI_STATUS_SIZE) :: status
+#endif
+
!----
! write seismograms in ASCII format
@@ -60,143 +72,172 @@
stop 'wrong component to save for seismograms'
endif
- do irec = 1,nrec
! only one seismogram if pressurs
- if(seismotype == 4) then
- number_of_components = 1
- else
- number_of_components = NDIM
- endif
+ if(seismotype == 4) then
+ number_of_components = 1
+ else
+ number_of_components = NDIM
+ endif
- do iorientation = 1,number_of_components
+ allocate(buffer_binary(NSTEP,number_of_components))
- if(iorientation == 1) then
- chn = 'BHX'
- else if(iorientation == 2) then
- chn = 'BHZ'
- else
- stop 'incorrect channel value'
- endif
+
+ if ( myrank == 0 ) then
+
+! delete the old files
+ open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown')
+ close(11,status='delete')
-! in case of pressure, use different abbreviation
- if(seismotype == 4) chn = 'PRE'
+ open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
+ close(11,status='delete')
+
+! write the new files
+ if(seismotype == 4) then
+ open(unit=12,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown',access='direct',recl=4*NSTEP)
+ else
+ open(unit=12,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown',access='direct',recl=4*NSTEP)
+ endif
-! create the name of the seismogram file for each slice
-! file name includes the name of the station, the network and the component
- length_station_name = len_trim(station_name(irec))
- length_network_name = len_trim(network_name(irec))
+ if(seismotype == 4) then
+ open(unit=13,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown',access='direct',recl=8*NSTEP)
+ else
+ open(unit=13,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8*NSTEP)
+ endif
-! check that length conforms to standard
- if(length_station_name < 1 .or. length_station_name > MAX_LENGTH_STATION_NAME) stop 'wrong length of station name'
+! no Z component seismogram if pressurs
+ if(seismotype /= 4) then
+ open(unit=14,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4*NSTEP)
+ open(unit=15,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8*NSTEP)
+
+ end if
+
+ end if
- if(length_network_name < 1 .or. length_network_name > MAX_LENGTH_NETWORK_NAME) stop 'wrong length of network name'
- write(sisname,"('OUTPUT_FILES/',a,'.',a,'.',a3,'.sem',a1)") station_name(irec)(1:length_station_name),&
- network_name(irec)(1:length_network_name),chn,component
-
-! save seismograms in text format with no subsampling.
-! Because we do not subsample the output, this can result in large files
-! if the simulation uses many time steps. However, subsampling the output
-! here would result in a loss of accuracy when one later convolves
-! the results with the source time function
- open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown')
-
-! make sure we never write more than the maximum number of time steps
-! subtract offset of the source to make sure travel time is correct
- do isample = 1,min(it,NSTEP)
- if(iorientation == 1) then
- write(11,*) sngl(dble(isample-1)*deltat - t0),' ',sngl(sisux(isample,irec))
+ irecloc = 0
+ do irec = 1,nrec
+
+ if ( myrank == 0 ) then
+
+ if ( which_proc_receiver(irec) == myrank ) then
+ irecloc = irecloc + 1
+ buffer_binary(:,1) = sisux(:,irecloc)
+ if ( number_of_components == 2 ) then
+ buffer_binary(:,2) = sisuz(:,irecloc)
+ end if
+
+#ifdef USE_MPI
else
- write(11,*) sngl(dble(isample-1)*deltat - t0),' ',sngl(sisuz(isample,irec))
- endif
- enddo
+ call MPI_RECV(buffer_binary(1,1),NSTEP,MPI_DOUBLE_PRECISION,&
+ which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
+ if ( number_of_components == 2 ) then
+ call MPI_RECV(buffer_binary(1,2),NSTEP,MPI_DOUBLE_PRECISION,&
+ which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
+ end if
+
+
+#endif
+ end if
+
+! write trace
+ do iorientation = 1,number_of_components
+
+ if(iorientation == 1) then
+ chn = 'BHX'
+ else if(iorientation == 2) then
+ chn = 'BHZ'
+ else
+ stop 'incorrect channel value'
+ endif
+
+ ! in case of pressure, use different abbreviation
+ if(seismotype == 4) chn = 'PRE'
+
+ ! create the name of the seismogram file for each slice
+ ! file name includes the name of the station, the network and the component
+ length_station_name = len_trim(station_name(irec))
+ length_network_name = len_trim(network_name(irec))
- close(11)
+ ! check that length conforms to standard
+ if(length_station_name < 1 .or. length_station_name > MAX_LENGTH_STATION_NAME) stop 'wrong length of station name'
+
+ if(length_network_name < 1 .or. length_network_name > MAX_LENGTH_NETWORK_NAME) stop 'wrong length of network name'
+
+ write(sisname,"('OUTPUT_FILES/',a,'.',a,'.',a3,'.sem',a1)") station_name(irec)(1:length_station_name),&
+ network_name(irec)(1:length_network_name),chn,component
+
+ ! save seismograms in text format with no subsampling.
+ ! Because we do not subsample the output, this can result in large files
+ ! if the simulation uses many time steps. However, subsampling the output
+ ! here would result in a loss of accuracy when one later convolves
+ ! the results with the source time function
+ open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown')
+
+ ! make sure we never write more than the maximum number of time steps
+ ! subtract offset of the source to make sure travel time is correct
+ do isample = 1,min(it,NSTEP)
+ if(iorientation == 1) then
+ write(11,*) sngl(dble(isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
+ else
+ write(11,*) sngl(dble(isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
+ endif
+ enddo
+
+ close(11)
+ end do
- enddo
+! write binary seismogram
+ write(12,rec=irec) sngl(buffer_binary(:,1))
+ write(13,rec=irec) sngl(buffer_binary(:,1))
+ if ( seismotype /= 4 ) then
+ write(14,rec=irec) sngl(buffer_binary(:,2))
+ write(15,rec=irec) sngl(buffer_binary(:,2))
+ end if
- enddo
+#ifdef USE_MPI
-!----
+ else
+ if ( which_proc_receiver(irec) == myrank ) then
+ irecloc = irecloc + 1
+ call MPI_SEND(sisux(1,irecloc),NSTEP,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ if ( number_of_components == 2 ) then
+ call MPI_SEND(sisuz(1,irecloc),NSTEP,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ end if
+ end if
-! write seismograms in single precision SEP binary format
+#endif
+
+ end if
-! X component
-
-! delete the old files
- open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown')
- close(11,status='delete')
-
- open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown')
- close(11,status='delete')
-
- open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown')
- close(11,status='delete')
-
- open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown')
- close(11,status='delete')
-
- irecord = 0
- do irec=1,nrec
- do isample=1,NSTEP
- irecord = irecord + 1
- buffer_binary_single(irecord) = sngl(sisux(isample,irec))
- buffer_binary_double(irecord) = sisux(isample,irec)
- enddo
enddo
-! write the new files
- if(seismotype == 4) then
- open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown',access='direct',recl=4*NSTEP*nrec)
- else
- open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown',access='direct',recl=4*NSTEP*nrec)
- endif
- write(11,rec=1) buffer_binary_single
- close(11)
+ close(12)
+ close(13)
+ if ( seismotype /= 4 ) then
+ close(14)
+ close(15)
+ end if
- if(seismotype == 4) then
- open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown',access='direct',recl=8*NSTEP*nrec)
- else
- open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8*NSTEP*nrec)
- endif
- write(11,rec=1) buffer_binary_double
- close(11)
+!----
-! Z component
-! delete the old files
- open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown')
- close(11,status='delete')
- open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
- close(11,status='delete')
-
-! no Z component seismogram if pressurs
- if(seismotype /= 4) then
-
- irecord = 0
- do irec=1,nrec
- do isample=1,NSTEP
- irecord = irecord + 1
- buffer_binary_single(irecord) = sngl(sisuz(isample,irec))
- buffer_binary_double(irecord) = sisuz(isample,irec)
- enddo
- enddo
-
-! write the new files
- open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4*NSTEP*nrec)
- write(11,rec=1) buffer_binary_single
- close(11)
-
- open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8*NSTEP*nrec)
- write(11,rec=1) buffer_binary_double
- close(11)
-
- endif
-
!----
-
+ if ( myrank == 0 ) then
+
! ligne de recepteurs pour Xsu
open(unit=11,file='OUTPUT_FILES/receiver_line_Xsu_XWindow',status='unknown')
@@ -242,6 +283,8 @@
write(11,*) '/bin/rm -f tempfile tempfile2'
close(11)
+end if
+
! formats
110 format('xwigb at xcur=',f8.2,'@n1=',i6,'@d1=',f15.8,'@f1=',f15.8,'@label1="Time@(s)"@label2="x@(m)"@n2=',i6,'@x2=')
More information about the cig-commits
mailing list