[cig-commits] r20711 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Sep 11 14:42:50 PDT 2012
Author: dkomati1
Date: 2012-09-11 14:42:49 -0700 (Tue, 11 Sep 2012)
New Revision: 20711
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
added MPI support to the calculation of total energy.
Also fixed a tiny bug that I introduced a few days ago in Postscript vector plots.
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90 2012-09-11 16:36:04 UTC (rev 20710)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90 2012-09-11 21:42:49 UTC (rev 20711)
@@ -48,13 +48,12 @@
xix,xiz,gammax,gammaz,jacobian,ibool, &
elastic,poroelastic,hprime_xx,hprime_zz, &
nspec,nglob_acoustic,nglob_elastic,nglob_poroelastic, &
- assign_external_model,it,deltat,t0,kmato,poroelastcoef,density, &
- porosity,tortuosity, &
+ assign_external_model,kmato,poroelastcoef,density,porosity,tortuosity, &
vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
anisotropic,anisotropy,wxgll,wzgll,numat, &
pressure_element,vector_field_element,e1,e11, &
potential_dot_acoustic,potential_dot_dot_acoustic, &
- ATTENUATION_VISCOELASTIC_SOLID,Mu_nu1,Mu_nu2,N_SLS,p_sv)
+ ATTENUATION_VISCOELASTIC_SOLID,Mu_nu1,Mu_nu2,N_SLS,p_sv,kinetic_energy,potential_energy)
! compute kinetic and potential energy in the solid (acoustic elements are excluded)
@@ -80,9 +79,6 @@
logical :: ATTENUATION_VISCOELASTIC_SOLID,p_sv
- integer :: it
- double precision :: t0,deltat
-
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
logical, dimension(nspec) :: elastic,poroelastic,anisotropic
@@ -395,9 +391,5 @@
enddo
- ! save kinetic, potential and total energy for this time step in external file
- write(IOUT_ENERGY,*) real(dble(it-1)*deltat - t0,4),real(kinetic_energy,4), &
- real(potential_energy,4),real(kinetic_energy + potential_energy,4)
-
end subroutine compute_energy
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90 2012-09-11 16:36:04 UTC (rev 20710)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90 2012-09-11 21:42:49 UTC (rev 20711)
@@ -2232,6 +2232,7 @@
endif
! draw the Stacey absorbing boundary line segment in different colors depending on its type
+ if ( myrank == 0 ) then
if(typeabs(inum) == IBOTTOM) then
write(24,*) '0 1 0 RG' ! Green
else if(typeabs(inum) == IRIGHT) then
@@ -2243,6 +2244,7 @@
else
call exit_MPI('Wrong typeabs() absorbing boundary code')
endif
+ endif
x1 = (coorg(1,knods(ideb,ispec))-xmin)*ratio_page + orig_x
z1 = (coorg(2,knods(ideb,ispec))-zmin)*ratio_page + orig_z
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-11 16:36:04 UTC (rev 20710)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-11 21:42:49 UTC (rev 20711)
@@ -362,6 +362,7 @@
include "constants.h"
#ifdef USE_MPI
include "mpif.h"
+ include "precision.h"
#endif
integer NSOURCES,i_source
@@ -550,6 +551,8 @@
double precision :: vpImin,vpImax,vpIImin,vpIImax
+ real(kind=CUSTOM_REAL) :: kinetic_energy,potential_energy,kinetic_energy_total,potential_energy_total
+
integer :: colors,numbers,subsamp_postscript,imagetype_postscript, &
NSTEP_BETWEEN_OUTPUT_INFO,seismotype,NSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP_BETWEEN_OUTPUT_IMAGES, &
NSTEP_BETWEEN_OUTPUT_WAVE_DUMPS,subsamp_seismos,imagetype_JPEG,imagetype_wavefield_dumps
@@ -4231,10 +4234,6 @@
endif
-#ifdef USE_MPI
- if(output_energy) stop 'error: energy calculation currently serial only, should add an MPI_REDUCE in parallel'
-#endif
-
!---- create a Gnuplot script to display the energy curve in log scale
if(output_energy .and. myrank == 0) then
close(IOUT_ENERGY)
@@ -7327,20 +7326,34 @@
endif
!---- compute kinetic and potential energy
- if(output_energy) &
+ if(output_energy) then
+
call compute_energy(displ_elastic,veloc_elastic, &
displs_poroelastic,velocs_poroelastic, &
displw_poroelastic,velocw_poroelastic, &
xix,xiz,gammax,gammaz,jacobian,ibool,elastic,poroelastic,hprime_xx,hprime_zz, &
nspec,nglob_acoustic,nglob_elastic,nglob_poroelastic, &
- assign_external_model,it,deltat,t0,kmato,poroelastcoef,density, &
- porosity,tortuosity, &
+ assign_external_model,kmato,poroelastcoef,density,porosity,tortuosity, &
vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
anisotropic,anisotropy,wxgll,wzgll,numat, &
pressure_element,vector_field_element,e1,e11, &
potential_dot_acoustic,potential_dot_dot_acoustic, &
- ATTENUATION_VISCOELASTIC_SOLID,Mu_nu1,Mu_nu2,N_SLS,p_sv)
+ ATTENUATION_VISCOELASTIC_SOLID,Mu_nu1,Mu_nu2,N_SLS,p_sv,kinetic_energy,potential_energy)
+#ifdef USE_MPI
+ call MPI_REDUCE(kinetic_energy, kinetic_energy_total, 1, CUSTOM_MPI_TYPE, MPI_SUM, 0, MPI_COMM_WORLD, ier)
+ call MPI_REDUCE(potential_energy, potential_energy_total, 1, CUSTOM_MPI_TYPE, MPI_SUM, 0, MPI_COMM_WORLD, ier)
+#else
+ kinetic_energy_total = kinetic_energy
+ potential_energy_total = potential_energy
+#endif
+
+! save kinetic, potential and total energy for this time step in external file
+ if(myrank == 0) write(IOUT_ENERGY,*) real(dble(it-1)*deltat - t0,4),real(kinetic_energy_total,4), &
+ real(potential_energy_total,4),real(kinetic_energy_total + potential_energy_total,4)
+
+ endif
+
!---- display time step and max of norm of displacement
if(mod(it,NSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
call check_stability(myrank,time,it,NSTEP,NOISE_TOMOGRAPHY, &
More information about the CIG-COMMITS
mailing list