[cig-commits] r8543 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:55:04 PST 2007
Author: walter
Date: 2007-12-07 15:55:03 -0800 (Fri, 07 Dec 2007)
New Revision: 8543
Added:
seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/plotpost.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
added calculation of elastic energy, taken from a very old version.
Also suppressed Nicolas' commented lines, starting with !!$
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:55:03 UTC (rev 8543)
@@ -26,7 +26,8 @@
LINK = $(F90)
-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
+#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
+LIB =
OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o
@@ -37,7 +38,7 @@
$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/construct_acoustic_surface.o $O/assemble_MPI.o
+ $O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_elastic_energy.o
default: clean meshfem2D specfem2D convolve_source_timefunction
@@ -130,6 +131,9 @@
$O/compute_gradient_attenuation.o: compute_gradient_attenuation.f90 constants.h
${F90} $(FLAGS_NOCHECK) -c -o $O/compute_gradient_attenuation.o compute_gradient_attenuation.f90
+$O/compute_elastic_energy.o: compute_elastic_energy.f90 constants.h
+ ${F90} $(FLAGS_NOCHECK) -c -o $O/compute_elastic_energy.o compute_elastic_energy.f90
+
$O/compute_vector_field.o: compute_vector_field.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/compute_vector_field.o compute_vector_field.f90
Modified: seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -19,13 +19,6 @@
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
@@ -129,33 +122,6 @@
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
Modified: seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -22,7 +22,6 @@
integer nspec
double precision xi_source,gamma_source
-!!!!!!!!!!!!!!!!! double precision Mxx,Myy,Mzz,Mxy,Mxz,Myz
double precision Mxx,Mzz,Mxz
double precision, dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
Added: seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -0,0 +1,152 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Dimitri Komatitsch
+! University of Pau, France
+!
+! (c) April 2007
+!
+!========================================================================
+
+ subroutine compute_elastic_energy(displ_elastic,veloc_elastic, &
+ xix,xiz,gammax,gammaz,jacobian,ibool,elastic,hprime_xx,hprime_zz, &
+ nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
+ vpext,vsext,rhoext,wxgll,wzgll,numat)
+
+! compute kinetic and potential energy in the solid (acoustic elements are excluded)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: nspec,npoin,numat
+
+ integer :: it
+ double precision :: t0,deltat
+
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+
+ logical, dimension(nspec) :: elastic
+
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
+
+ integer, dimension(nspec) :: kmato
+
+ logical :: assign_external_model
+
+ double precision, dimension(numat) :: density
+ double precision, dimension(4,numat) :: elastcoef
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
+
+ double precision, dimension(NDIM,npoin) :: displ_elastic,veloc_elastic
+
+! Gauss-Lobatto-Legendre points and weights
+ double precision, dimension(NGLLX) :: wxgll
+ double precision, dimension(NGLLZ) :: wzgll
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLLX,NGLLX) :: hprime_xx
+ double precision, dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! local variables
+ integer :: i,j,k,ispec
+
+! spatial derivatives
+ double precision :: dux_dxi,dux_dgamma,duz_dxi,duz_dgamma
+ double precision :: dux_dxl,duz_dxl,dux_dzl,duz_dzl
+
+! jacobian
+ double precision :: xixl,xizl,gammaxl,gammazl,jacobianl
+
+ double precision :: kinetic_energy,potential_energy
+ double precision :: cpl,csl,rhol,mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed
+
+ kinetic_energy = ZERO
+ potential_energy = ZERO
+
+! loop over spectral elements
+ do ispec = 1,nspec
+
+!---
+!--- elastic spectral element
+!---
+ if(elastic(ispec)) then
+
+! get relaxed elastic parameters of current spectral element
+ lambdal_relaxed = elastcoef(1,kmato(ispec))
+ mul_relaxed = elastcoef(2,kmato(ispec))
+ lambdalplus2mul_relaxed = elastcoef(3,kmato(ispec))
+ rhol = density(kmato(ispec))
+
+! first double loop over GLL points to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+!--- if external medium, get elastic parameters of current grid point
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ mul_relaxed = rhol*csl*csl
+ lambdal_relaxed = rhol*cpl*cpl - TWO*mul_relaxed
+ lambdalplus2mul_relaxed = lambdal_relaxed + TWO*mul_relaxed
+ endif
+
+! derivative along x and along z
+ dux_dxi = ZERO
+ duz_dxi = ZERO
+
+ dux_dgamma = ZERO
+ duz_dgamma = ZERO
+
+! first double loop over GLL points to compute and store gradients
+! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(k,i)
+ duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(k,i)
+ dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(k,j)
+ duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(k,j)
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+ jacobianl = jacobian(i,j,ispec)
+
+! derivatives of displacement
+ dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+
+ duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
+ duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
+
+! compute potential energy
+ potential_energy = potential_energy + (lambdalplus2mul_relaxed*dux_dxl**2 &
+ + lambdalplus2mul_relaxed*duz_dzl**2 &
+ + two*lambdal_relaxed*dux_dxl*duz_dzl + mul_relaxed*(dux_dzl + duz_dxl)**2)*wxgll(i)*wzgll(j)*jacobianl
+
+! compute kinetic energy
+ kinetic_energy = kinetic_energy + &
+ rhol*(veloc_elastic(1,ibool(i,j,ispec))**2 + veloc_elastic(2,ibool(i,j,ispec))**2) *wxgll(i)*wzgll(j)*jacobianl
+
+ enddo
+ enddo
+
+ endif
+
+ enddo
+
+! do not forget to divide by two at the end
+ kinetic_energy = kinetic_energy / TWO
+ potential_energy = potential_energy / TWO
+
+! save kinetic, potential and total energy for this time step in external file
+ write(IENERGY,*) sngl(dble(it-1)*deltat - t0),sngl(kinetic_energy), &
+ sngl(potential_energy),sngl(kinetic_energy + potential_energy)
+
+ end subroutine compute_elastic_energy
+
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:55:03 UTC (rev 8543)
@@ -4,6 +4,11 @@
! the code does NOT work if NGLLZ /= NGLLX because it then cannot handle a non-structured mesh
integer, parameter :: NGLLZ = NGLLX
+! compute and output elastic energy (slows down the code significantly)
+ logical, parameter :: OUTPUT_ELASTIC_ENERGY = .false.
+! output file for energy
+ integer, parameter :: IENERGY = 77
+
! select fast (Paul Fischer) or slow (topology only) global numbering algorithm
logical, parameter :: FAST_NUMBERING = .true.
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -50,10 +50,6 @@
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
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -347,7 +347,6 @@
end do
end if
- !nnodes = (nx+1) * (nz+1)
end if
! read absorbing boundaries parameters
@@ -873,20 +872,12 @@
!!$ print *,'Grid saved in Gnuplot format...'
!!$ print *
-
!*****************************
! Partitionning
!*****************************
allocate(part(0:nelmnts-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, we work on a subarray of elmnts, which represents the elements with for nodes only
! construction of the graph
if ( ngnod == 9 ) then
@@ -1066,9 +1057,6 @@
write(15,*) nb_materials,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
- !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, &
@@ -1122,9 +1110,6 @@
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, &
Modified: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -468,8 +468,6 @@
allocate(tab_size_interfaces(0:ninterfaces))
tab_size_interfaces(:) = 0
-! print *, 'num_interface', ninterfaces
-
num_interface = 0
num_edge = 0
@@ -489,14 +487,12 @@
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
@@ -526,7 +522,6 @@
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
@@ -542,7 +537,7 @@
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
+ print *, "Error while building interfaces!", ncommon_nodes
end if
num_edge = num_edge + 1
end if
@@ -930,7 +925,6 @@
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)
@@ -1018,133 +1012,6 @@
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
@@ -1375,7 +1242,6 @@
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'
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.F90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -2197,9 +2197,6 @@
!
!--- draw free surface with a thick color line
!
-#ifdef USE_MPI
-! call MPI_ALLREDUCE(anyabs, anyabs_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
-#endif
if ( myrank == 0 ) then
write(24,*) '%'
@@ -2223,26 +2220,6 @@
do inum = 1,nelem_acoustic_surface
ispec = acoustic_edges(1,inum)
-!!$ do iedge = 1,4
-!!$
-!!$ if(codeabs(iedge,inum) /= 0) then
-!!$
-!!$ if(iedge == ITOP) then
-!!$ ideb = 3
-!!$ ifin = 4
-!!$ else if(iedge == IBOTTOM) then
-!!$ ideb = 1
-!!$ ifin = 2
-!!$ else if(iedge == ILEFT) then
-!!$ ideb = 4
-!!$ ifin = 1
-!!$ else if(iedge == IRIGHT) then
-!!$ ideb = 2
-!!$ ifin = 3
-!!$ else
-!!$ call exit_MPI('Wrong absorbing boundary code')
-!!$ endif
-
x1 = (coorg(1,acoustic_edges(3,inum))-xmin)*ratio_page + orig_x
z1 = (coorg(2,acoustic_edges(3,inum))-zmin)*ratio_page + orig_z
x2 = (coorg(1,acoustic_edges(4,inum))-xmin)*ratio_page + orig_x
@@ -2261,9 +2238,6 @@
coorg_send(4,buffer_offset) = z2
end if
-!!$ endif
-!!$ enddo
-
enddo
end if
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -244,7 +244,7 @@
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
@@ -259,7 +259,7 @@
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
@@ -269,7 +269,7 @@
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 :: max_ibool_interfaces_size_acoustic, max_ibool_interfaces_size_elastic
integer, dimension(:,:), allocatable :: acoustic_surface
integer, dimension(:,:), allocatable :: acoustic_edges
@@ -282,7 +282,7 @@
integer :: nrecloc, irecloc
integer, dimension(:), allocatable :: recloc, which_proc_receiver
-
+
character(len=256) :: filename
!***********************************************************************
@@ -295,19 +295,19 @@
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')
@@ -472,16 +472,6 @@
allocate(jbegin_right(nelemabs))
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
-
!
!---- print element group main parameters
!
@@ -556,7 +546,7 @@
!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))
@@ -577,10 +567,10 @@
end do
end do
print *, 'end read the interfaces', myrank
-
+
end if
-
-
+
+
!
!---- read absorbing boundary data
!
@@ -588,7 +578,7 @@
if(anyabs) then
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)
+ jbegin_right(inum), jend_right(inum), ibegin_top(inum), iend_top(inum), jbegin_left(inum), jend_left(inum)
if(numabsread < 1 .or. numabsread > nspec) call exit_MPI('Wrong absorbing element number')
numabs(inum) = numabsread
codeabs(IBOTTOM,inum) = codeabsread(1)
@@ -613,12 +603,6 @@
allocate(acoustic_surface(5,nelem_acoustic_surface))
call construct_acoustic_surface ( nspec, ngnod, knods, nelem_acoustic_surface, &
acoustic_edges, acoustic_surface)
-!!$ read(IIN,*) numacoustread,iedgeacoustread
-!!$ if(numacoustread < 1 .or. numacoustread > nspec) call eixt_MPI('Wrong acoustic free surface element number')
-!!$ if(iedgeacoustread < 1 .or. iedgeacoustread > NEDGES) call exit_MPI('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
@@ -640,7 +624,7 @@
allocate(fluid_solid_acoustic_iedge(1))
allocate(fluid_solid_elastic_ispec(1))
allocate(fluid_solid_elastic_iedge(1))
-
+
end if
@@ -830,7 +814,7 @@
! for attenuation
if(TURN_ANISOTROPY_ON .and. TURN_ATTENUATION_ON) then
call exit_MPI('cannot have anisotropy and attenuation both turned on in current version')
- end if
+ end if
!
!---- define coefficients of the Newmark time scheme
!
@@ -860,7 +844,7 @@
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, &
@@ -994,7 +978,7 @@
if ( nproc > 1 ) then
! preparing for MPI communications
call prepare_assemble_MPI (myrank,nspec,ibool, &
- knods, ngnod, &
+ knods, ngnod, &
npoin, elastic, &
ninterface, max_interface_size, &
my_neighbours, my_nelmnts_neighbours, my_interfaces, &
@@ -1024,7 +1008,7 @@
tab_requests_send_recv_acoustic, &
inum_interfaces_acoustic &
)
-
+
! creating mpi non-blocking persistent communications for elastic elements
call create_MPI_requests_SEND_RECV_elastic(myrank, &
ninterface, ninterface_elastic, &
@@ -1078,7 +1062,7 @@
zmin_color_image_loc = minval(coord(2,:))
zmax_color_image_loc = maxval(coord(2,:))
-! global values
+! global values
xmin_color_image = xmin_color_image_loc
xmax_color_image = xmax_color_image_loc
zmin_color_image = zmin_color_image_loc
@@ -1102,7 +1086,7 @@
! 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))
-
+
! convert pixel sizes to even numbers because easier to reduce size, create MPEG movies in postprocessing
NX_IMAGE_color = 2 * (NX_IMAGE_color / 2)
NZ_IMAGE_color = 2 * (NZ_IMAGE_color / 2)
@@ -1110,10 +1094,10 @@
! check that image size is not too big
if ( NX_IMAGE_color > 99999 ) then
call exit_MPI('output image too big : NX_IMAGE_color > 99999.')
- end if
+ end if
if ( NZ_IMAGE_color > 99999 ) then
call exit_MPI('output image too big : NZ_IMAGE_color > 99999.')
- end if
+ end if
! allocate an array for image data
allocate(image_color_data(NX_IMAGE_color,NZ_IMAGE_color))
@@ -1134,24 +1118,24 @@
! cheking which pixels are inside each elements
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))
end do
-
+
! avoid working on the whole pixel grid
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
-
+
! checking if the pixel is inside the element (must be a convex quadrilateral)
call is_in_convex_quadrilateral( elmnt_coords, i_coord, j_coord, pixel_is_in)
@@ -1165,46 +1149,46 @@
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
call exit_MPI('Error in detecting pixel for color image')
-
+
end if
nb_pixel_loc = nb_pixel_loc + 1
-
+
end if
-
+
end do
end do
end do
! creating and filling array num_pixel_loc with the positions of each colored pixel owned by the local process (useful for parallel jobs)
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
-
+
! filling array iglob_image_color, containing info on which process owns which pixels.
#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)
-
+
+ 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(:))))
@@ -1213,7 +1197,7 @@
allocate(data_pixel_send(nb_pixel_loc))
if ( nproc > 1 ) then
if ( myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV(num_pixel_recv(1,iproc+1),nb_pixel_per_proc(iproc+1), MPI_INTEGER, &
@@ -1222,13 +1206,13 @@
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
end if
@@ -1337,11 +1321,11 @@
write(55,*) sngl(time),sngl(source_time_function(it)),sngl(time-t0)
endif
enddo
-
+
if ( myrank == 0 ) then
close(55)
endif
-
+
source_time_function(:) = source_time_function(:) / nb_proc_source
else
@@ -1355,7 +1339,7 @@
coupled_acoustic_elastic = any_acoustic .and. any_elastic
! fluid/solid edge detection
-! the two elements (fluid and solid) forming an edge are already known (computed in meshfem2D),
+! the two elements (fluid and solid) forming an edge are already known (computed in meshfem2D),
! the common nodes forming the edge are computed here
if(coupled_acoustic_elastic) then
print *
@@ -1438,8 +1422,8 @@
endif
enddo
-
+
!if(num_fluid_solid_edges /= num_fluid_solid_edges_alloc) call exit_MPI('error in creation of arrays for fluid/solid matching')
! make sure fluid/solid matching has been perfectly detected: check that the grid points
@@ -1487,19 +1471,6 @@
endif
-
-! default values for acoustic absorbing edges
-!!$ ibegin_bottom(:) = 1
-!!$ ibegin_top(:) = 1
-!!$
-!!$ 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
@@ -1520,7 +1491,7 @@
! if acoustic absorbing element and acoustic/elastic coupled element is the same
if(ispec_acoustic == ispec) then
-
+
if(iedge_acoustic == IBOTTOM) then
jbegin_left(ispecabs) = 2
jbegin_right(ispecabs) = 2
@@ -1549,6 +1520,12 @@
endif
+#ifdef USE_MPI
+ if(OUTPUT_ELASTIC_ENERGY) stop 'energy calculation only serial right now, should add an MPI_REDUCE in parallel'
+#endif
+! open the file in which we will store the energy curve
+ if(OUTPUT_ELASTIC_ENERGY) open(unit=IENERGY,file='energy.gnu',status='unknown')
+
!
!---- s t a r t t i m e i t e r a t i o n s
!
@@ -1575,7 +1552,7 @@
! *********************************************************
do it = 1,NSTEP
-
+
! compute current time
time = (it-1)*deltat
@@ -1699,7 +1676,7 @@
tab_requests_send_recv_acoustic, &
buffer_send_faces_vector_acoustic, &
buffer_recv_faces_vector_acoustic &
- )
+ )
end if
#endif
@@ -1718,7 +1695,7 @@
potential_acoustic,acoustic_surface, &
ibool,nelem_acoustic_surface,npoin,nspec)
endif
-
+
! *********************************************************
! ************* main solver for the elastic elements
! *********************************************************
@@ -1825,7 +1802,7 @@
tab_requests_send_recv_elastic, &
buffer_send_faces_vector_elastic, &
buffer_recv_faces_vector_elastic &
- )
+ )
end if
#endif
@@ -1840,6 +1817,13 @@
veloc_elastic = veloc_elastic + deltatover2*accel_elastic
endif
+!---- compute kinetic and potential energy
+!
+ if(OUTPUT_ELASTIC_ENERGY) &
+ call compute_elastic_energy(displ_elastic,veloc_elastic, &
+ xix,xiz,gammax,gammaz,jacobian,ibool,elastic,hprime_xx,hprime_zz, &
+ nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
+ vpext,vsext,rhoext,wxgll,wzgll,numat)
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
@@ -1870,9 +1854,9 @@
! loop on all the receivers to compute and store the seismograms
do irecloc = 1,nrecloc
-
+
irec = recloc(irecloc)
-
+
ispec = ispec_selected_rec(irec)
! compute pressure in this element if needed
@@ -2078,64 +2062,64 @@
endif
image_color_data(:,:) = 0.d0
-
+
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
-
+
! assembling array image_color_data on process zero for color output
-#ifdef USE_MPI
+#ifdef USE_MPI
if ( nproc > 1 ) then
if ( myrank == 0 ) then
-
+
do iproc = 1, nproc-1
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)
-#endif
+#endif
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
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
@@ -2162,6 +2146,18 @@
enddo ! end of the main time loop
+!---- close energy file and create a gnuplot script to display it
+ if(OUTPUT_ELASTIC_ENERGY) then
+ close(IENERGY)
+ open(unit=IENERGY,file='plotenergy',status='unknown')
+ write(IENERGY,*) 'set term postscript landscape color solid "Helvetica" 22'
+ write(IENERGY,*) 'set output "energy.ps"'
+ write(IENERGY,*) 'set xlabel "Time (s)"'
+ write(IENERGY,*) 'set ylabel "Energy (J)"'
+ write(IENERGY,*) 'plot "energy.gnu" us 1:4 t ''Total Energy'' w l 1, "energy.gnu" us 1:3 t ''Potential Energy'' w l 2'
+ close(IENERGY)
+ endif
+
! print exit banner
call datim(simulation_title)
@@ -2244,17 +2240,17 @@
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
+ 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)
@@ -2271,11 +2267,10 @@
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-07-02 21:50:08 UTC (rev 8542)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-12-07 23:55:03 UTC (rev 8543)
@@ -23,13 +23,13 @@
include "mpif.h"
#endif
- integer nrec,NSTEP,it,seismotype
- double precision t0,deltat
+ integer :: nrec,NSTEP,it,seismotype
+ double precision :: t0,deltat
- integer, intent(in) :: nrecloc,myrank
- integer, dimension(nrec),intent(in) :: which_proc_receiver
+ integer, intent(in) :: nrecloc,myrank
+ integer, dimension(nrec),intent(in) :: which_proc_receiver
- double precision, dimension(NSTEP,nrecloc), intent(in) :: sisux,sisuz
+ double precision, dimension(NSTEP,nrecloc), intent(in) :: sisux,sisuz
double precision st_xval(nrec)
More information about the cig-commits
mailing list