[cig-commits] r8555 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:56:06 PST 2007
Author: walter
Date: 2007-12-07 15:56:05 -0800 (Fri, 07 Dec 2007)
New Revision: 8555
Added:
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
Removed:
seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
added calculation of kinetic and potential energy in the fluid (we had it in the solid only before)
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-07-04 20:20:00 UTC (rev 8554)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:56:05 UTC (rev 8555)
@@ -38,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/compute_elastic_energy.o
+ $O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_energy.o
default: clean meshfem2D specfem2D convolve_source_timefunction
@@ -131,8 +131,8 @@
$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_energy.o: compute_energy.f90 constants.h
+ ${F90} $(FLAGS_NOCHECK) -c -o $O/compute_energy.o compute_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
Deleted: seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90 2007-07-04 20:20:00 UTC (rev 8554)
+++ seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90 2007-12-07 23:56:05 UTC (rev 8555)
@@ -1,152 +0,0 @@
-
-!========================================================================
-!
-! 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(i,k)
- duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
- dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
- 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
-
Copied: seismo/2D/SPECFEM2D/trunk/compute_energy.f90 (from rev 8554, seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90)
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_elastic_energy.f90 2007-07-04 20:20:00 UTC (rev 8554)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-12-07 23:56:05 UTC (rev 8555)
@@ -0,0 +1,217 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Dimitri Komatitsch
+! University of Pau, France
+!
+! (c) April 2007
+!
+!========================================================================
+
+ subroutine compute_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, &
+ pressure_element,vector_field_element,e1_mech1,e11_mech1,e1_mech2,e11_mech2, &
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+
+! compute kinetic and potential energy in the solid (acoustic elements are excluded)
+
+ implicit none
+
+ include "constants.h"
+
+! vector field in an element
+ double precision, dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
+
+! pressure in an element
+ double precision, dimension(NGLLX,NGLLX) :: pressure_element
+
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: e1_mech1,e11_mech1,e1_mech2,e11_mech2
+
+ double precision, dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
+
+ logical :: TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+
+ 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,kappal
+
+ 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))
+
+! double loop over GLL points
+ 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(i,k)
+ duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ 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 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 / TWO
+
+! 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 / TWO
+
+ enddo
+ enddo
+
+!---
+!--- acoustic spectral element
+!---
+ else
+
+! for the definition of potential energy in an acoustic fluid, see for instance
+! equation (23) of M. Maess et al., Journal of Sound and Vibration 296 (2006) 264-276
+
+! in case of an acoustic medium, a displacement potential Chi is used as in Chaljub and Valette,
+! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
+! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
+! This permits acoustic-elastic coupling based on a non-iterative time scheme.
+! Displacement is then: u = grad(Chi)
+! Velocity is then: v = grad(Chi_dot) (Chi_dot being the time derivative of Chi)
+! and pressure is: p = - rho * Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
+
+! compute pressure in this element
+ call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
+ numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1_mech1,e11_mech1, &
+ e1_mech2,e11_mech2,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
+
+! compute velocity vector field in this element
+ call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+
+! get density of current spectral element
+ lambdal_relaxed = elastcoef(1,kmato(ispec))
+ mul_relaxed = elastcoef(2,kmato(ispec))
+ rhol = density(kmato(ispec))
+ kappal = lambdal_relaxed + TWO*mul_relaxed/3.d0
+ cpl = sqrt((kappal + 4.d0*mul_relaxed/3.d0)/rhol)
+
+! double loop over GLL points
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+!--- if external medium, get density of current grid point
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+
+! compute kinetic energy
+ kinetic_energy = kinetic_energy + &
+ rhol*(vector_field_element(1,i,j)**2 + &
+ vector_field_element(2,i,j)**2) *wxgll(i)*wzgll(j)*jacobianl / TWO
+
+! compute potential energy
+ potential_energy = potential_energy + (pressure_element(i,j)**2)*wxgll(i)*wzgll(j)*jacobianl / (TWO * rhol * cpl**2)
+
+ enddo
+ enddo
+
+ endif
+
+ enddo
+
+! 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_energy
+
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-07-04 20:20:00 UTC (rev 8554)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:56:05 UTC (rev 8555)
@@ -4,8 +4,8 @@
! 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.
+! compute and output acoustic and elastic energy (slows down the code significantly)
+ logical, parameter :: OUTPUT_ENERGY = .false.
! output file for energy
integer, parameter :: IENERGY = 77
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-07-04 20:20:00 UTC (rev 8554)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:56:05 UTC (rev 8555)
@@ -110,8 +110,10 @@
! for seismograms
double precision, dimension(:,:), allocatable :: sisux,sisuz
+
! vector field in an element
double precision, dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
+
! pressure in an element
double precision, dimension(NGLLX,NGLLX) :: pressure_element
@@ -1197,7 +1199,7 @@
end do
end do
-! creating and filling array num_pixel_loc with the positions of each colored
+! 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))
@@ -1560,10 +1562,10 @@
endif
#ifdef USE_MPI
- if(OUTPUT_ELASTIC_ENERGY) stop 'energy calculation only serial right now, should add an MPI_REDUCE in parallel'
+ if(OUTPUT_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')
+ if(OUTPUT_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
@@ -1874,11 +1876,13 @@
!---- compute kinetic and potential energy
!
- if(OUTPUT_ELASTIC_ENERGY) &
- call compute_elastic_energy(displ_elastic,veloc_elastic, &
+ if(OUTPUT_ENERGY) &
+ call compute_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)
+ vpext,vsext,rhoext,wxgll,wzgll,numat, &
+ pressure_element,vector_field_element,e1_mech1,e11_mech1,e1_mech2,e11_mech2, &
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON)
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
@@ -2197,7 +2201,7 @@
enddo ! end of the main time loop
!---- close energy file and create a gnuplot script to display it
- if(OUTPUT_ELASTIC_ENERGY) then
+ if(OUTPUT_ENERGY) then
close(IENERGY)
open(unit=IENERGY,file='plotenergy',status='unknown')
write(IENERGY,*) 'set term postscript landscape color solid "Helvetica" 22'
More information about the cig-commits
mailing list