[cig-commits] r19276 - seismo/1D/SPECFEM1D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Dec 6 16:29:58 PST 2011
Author: dkomati1
Date: 2011-12-06 16:29:58 -0800 (Tue, 06 Dec 2011)
New Revision: 19276
Added:
seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu
Modified:
seismo/1D/SPECFEM1D/trunk/diffusion.f90
seismo/1D/SPECFEM1D/trunk/wave.f90
Log:
renamed some variables for clarity, added a few comments, and added a Gnuplot script to display the solution
Modified: seismo/1D/SPECFEM1D/trunk/diffusion.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/diffusion.f90 2011-12-07 00:05:15 UTC (rev 19275)
+++ seismo/1D/SPECFEM1D/trunk/diffusion.f90 2011-12-07 00:29:58 UTC (rev 19276)
@@ -66,7 +66,7 @@
double precision, parameter :: THERMALCONDUCTIVITY = 10.0d-01 ! cal/m/s/K
double precision, parameter :: HEATCAPACITY = 0.3d+03 ! cal/kg/K
- integer ispec,i,j,iglob,itime
+ integer ispec,i,j,iglob,it
! Gauss-Lobatto-Legendre points of integration
double precision, dimension(NGLL) :: xigll
@@ -87,7 +87,7 @@
double precision, dimension(NGLL,NSPEC) :: rho,heat_capacity,thermal_conductivity
! Jacobian `matrix' and Jacobian
- double precision, dimension(NGLL,NSPEC) :: dxidx,jacobian
+ double precision, dimension(NGLL,NSPEC) :: dxi_dx,jacobian
! local mass matrix
double precision mass_local
@@ -112,7 +112,7 @@
double precision temperature_1,temperature_NGLOB
! derivatives
- double precision dtdx,flux,templ,temp(NGLL)
+ double precision dt_dx,flux,templ,temp(NGLL)
! movie
character(len=50) moviefile
@@ -133,7 +133,7 @@
rho(i,ispec) = DENSITY
thermal_conductivity(i,ispec) = THERMALCONDUCTIVITY
heat_capacity(i,ispec) = HEATCAPACITY
- dxidx(i,ispec) = 2. / (x2(ispec)-x1(ispec))
+ dxi_dx(i,ispec) = 2. / (x2(ispec)-x1(ispec)) ! this is d(xi) / dx
jacobian(i,ispec) = (x2(ispec)-x1(ispec)) / 2.
enddo
enddo
@@ -147,7 +147,7 @@
enddo
enddo
-! get the global grid points
+! compute the position of the global grid points
do ispec = 1,NSPEC
do i = 1,NGLL
iglob = ibool(i,ispec)
@@ -155,7 +155,7 @@
enddo
enddo
-! calculate the global mass matrix
+! calculate the assembled global mass matrix
mass_global(:) = 0.
do ispec = 1,NSPEC
do i = 1,NGLL
@@ -169,7 +169,7 @@
dh = LENGTH/dble(NGLOB-1)
diffusivity = THERMALCONDUCTIVITY/(HEATCAPACITY*DENSITY)
time_step = 0.5*dh*dh/diffusivity
-! print *,'time step estimate: ',time_step,' seconds'
+! print *,'time step estimate: ',time_step,' seconds'
if(FIXED_BC) then
! set up the temperatures at the ends
@@ -190,7 +190,7 @@
temperature(:) = 0.
grad_temperature(:) = 0.
- do itime = 1,NSTEP
+ do it = 1,NSTEP
! update temperature
temperature(:) = temperature(:) + deltatover2*grad_temperature(:)
@@ -199,18 +199,18 @@
do ispec = 1,NSPEC
do i = 1,NGLL
-! get dtdx
+! compute dt_dx
templ = 0.
do j = 1,NGLL
iglob = ibool(j,ispec)
templ = templ + temperature(iglob)*hprime(i,j)
enddo
- dtdx = templ*dxidx(i,ispec)
+ dt_dx = templ*dxi_dx(i,ispec)
! heat flux
- flux = -thermal_conductivity(i,ispec)*dtdx
+ flux = -thermal_conductivity(i,ispec)*dt_dx
- temp(i) = jacobian(i,ispec)*flux*dxidx(i,ispec)
+ temp(i) = jacobian(i,ispec)*flux*dxi_dx(i,ispec)
enddo ! first loop over the GLL points
do i = 1,NGLL
@@ -235,8 +235,8 @@
iglob = ibool(i,1)
templ = templ + temperature(iglob)*hprime(1,i)
enddo
- dtdx = templ*dxidx(1,1)
- flux_1 = -thermal_conductivity(1,1)*dtdx
+ dt_dx = templ*dxi_dx(1,1)
+ flux_1 = -thermal_conductivity(1,1)*dt_dx
! right side
temperature(NGLOB) = temperature_NGLOB
templ = 0.
@@ -244,8 +244,8 @@
iglob = ibool(i,NSPEC)
templ = templ + temperature(iglob)*hprime(NGLL,i)
enddo
- dtdx = templ*dxidx(NGLL,NSPEC)
- flux_NGLOB = -thermal_conductivity(NGLL,NSPEC)*dtdx
+ dt_dx = templ*dxi_dx(NGLL,NSPEC)
+ flux_NGLOB = -thermal_conductivity(NGLL,NSPEC)*dt_dx
endif
! add in the end fluxes
@@ -259,8 +259,8 @@
temperature(:) = temperature(:) + deltatover2*grad_temperature(:)
! write out snapshots
- if(mod(itime-1,1000) == 0) then
- write(moviefile,"('snapshot',i5.5)") itime
+ if(mod(it-1,1000) == 0) then
+ write(moviefile,"('snapshot',i5.5)") it
open(unit=10,file=moviefile,status='unknown')
do iglob = 1,NGLOB
write(10,*) sngl(x(iglob)),sngl(temperature(iglob))
@@ -271,3 +271,4 @@
enddo ! end time loop
end program diffusion
+
Added: seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu
===================================================================
--- seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu (rev 0)
+++ seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu 2011-12-07 00:29:58 UTC (rev 19276)
@@ -0,0 +1,37 @@
+# Gnuplot script to plot all the snapshots
+
+set term x11
+
+set xrange [0:3000]
+set yrange [-1.5:+1.5]
+
+plot "snapshot00001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot01001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot02001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot03001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot04001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot05001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot06001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot07001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot08001" w l lc 1
+pause -1 "hit any key..."
+
+plot "snapshot09001" w l lc 1
+pause -1 "hit any key..."
+
Modified: seismo/1D/SPECFEM1D/trunk/wave.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/wave.f90 2011-12-07 00:05:15 UTC (rev 19275)
+++ seismo/1D/SPECFEM1D/trunk/wave.f90 2011-12-07 00:29:58 UTC (rev 19276)
@@ -68,7 +68,7 @@
! pi
double precision, parameter :: PI = 3.141592653589793d0
- integer ispec,i,j,iglob,itime
+ integer ispec,i,j,iglob,it
! Gauss-Lobatto-Legendre points of integration
double precision, dimension(NGLL) :: xigll
@@ -89,7 +89,7 @@
double precision, dimension(NGLL,NSPEC) :: rho,mu
! Jacobian `matrix' and Jacobian
- double precision, dimension(NGLL,NSPEC) :: dxidx,jacobian
+ double precision, dimension(NGLL,NSPEC) :: dxi_dx,jacobian
! local mass matrix
double precision mass_local
@@ -104,7 +104,7 @@
integer, dimension(NGLL,NSPEC) :: ibool
! time marching
- double precision dh,v,courant,time_step
+ double precision dh,v,courant_CFL,time_step
double precision deltat,deltatover2,deltatsqover2
! source
@@ -117,7 +117,7 @@
double precision seismogram(NSTEP)
! derivatives
- double precision dudx,sigma,templ,temp(NGLL)
+ double precision du_dxi,epsilon,sigma,templ,temp(NGLL)
! movie
character(len=50) moviefile
@@ -137,7 +137,7 @@
do i = 1,NGLL
rho(i,ispec) = DENSITY
mu(i,ispec) = RIGIDITY
- dxidx(i,ispec) = 2. / (x2(ispec)-x1(ispec))
+ dxi_dx(i,ispec) = 2. / (x2(ispec)-x1(ispec)) ! this is d(xi) / dx
jacobian(i,ispec) = (x2(ispec)-x1(ispec)) / 2.
enddo
enddo
@@ -151,7 +151,7 @@
enddo
enddo
-! get the global grid points
+! compute the position of the global grid points
do ispec = 1,NSPEC
do i = 1,NGLL
iglob = ibool(i,ispec)
@@ -159,7 +159,7 @@
enddo
enddo
-! calculate the global mass matrix
+! calculate the assembled global mass matrix
mass_global(:) = 0.
do ispec = 1,NSPEC
do i = 1,NGLL
@@ -172,9 +172,9 @@
! estimate the time step
dh = LENGTH/dble(NGLOB-1)
v = dsqrt(RIGIDITY/DENSITY)
- courant = 0.2
- time_step = courant*dh/v
-! print *,'time step estimate: ',time_step,' seconds'
+ courant_CFL = 0.2
+ time_step = courant_CFL*dh/v
+! print *,'time step estimate: ',time_step,' seconds'
! set the source
ispec_source = (NSPEC+1)/2
@@ -199,9 +199,9 @@
displ(iglob) = dsin(PI*x(iglob)/LENGTH)
enddo
- do itime = 1,NSTEP
+ do it = 1,NSTEP
-! `predictor' update displacement using finite-difference time scheme (Newmark)
+! `predictor' update displacement using explicit finite-difference time scheme (Newmark)
displ(:) = displ(:) + deltat*veloc(:) + deltatsqover2*accel(:)
veloc(:) = veloc(:) + deltatover2*accel(:)
accel(:) = 0.
@@ -209,18 +209,20 @@
do ispec = 1,NSPEC
do i = 1,NGLL
-! get dudx
- templ = 0.
+! compute d(u) / d(xi)
+ du_dxi = 0.
do j = 1,NGLL
iglob = ibool(j,ispec)
- templ = templ + displ(iglob)*hprime(i,j)
+ du_dxi = du_dxi + displ(iglob)*hprime(i,j)
enddo
- dudx = templ*dxidx(i,ispec)
+! strain
+ epsilon = du_dxi*dxi_dx(i,ispec)
+
! stress
- sigma = mu(i,ispec)*dudx
+ sigma = mu(i,ispec)*epsilon
- temp(i) = jacobian(i,ispec)*sigma*dxidx(i,ispec)
+ temp(i) = jacobian(i,ispec)*sigma*dxi_dx(i,ispec)
enddo ! first loop over the GLL points
do i = 1,NGLL
@@ -229,17 +231,19 @@
templ = templ + temp(j)*hprime(j,i)*wgll(j)
enddo
-! `corrector' update acceleration
+! `corrector' update of acceleration in the Newmark scheme
+! the minus sign comes from the integration by part done in the weak formulation of the equations
iglob = ibool(i,ispec)
accel(iglob) = accel(iglob) - templ
+
enddo ! second loop over the GLL points
enddo ! end loop over all spectral elements
! add source at global level
-! iglob_source = ibool(ispec_source,i_source)
-! stf = source_time_function(dble(itime-1)*DT-hdur,hdur)
-! accel(iglob_source) = accel(iglob_source) + stf*source_amp
+! iglob_source = ibool(ispec_source,i_source)
+! stf = source_time_function(dble(it-1)*DT-hdur,hdur)
+! accel(iglob_source) = accel(iglob_source) + stf*source_amp
! fixed boundary conditions at global level
if(FIXED_BC) then
@@ -254,8 +258,8 @@
veloc(:) = veloc(:) + deltatover2*accel(:)
! write out snapshots
- if(mod(itime-1,1000) == 0) then
- write(moviefile,"('snapshot',i5.5)") itime
+ if(mod(it-1,1000) == 0) then
+ write(moviefile,"('snapshot',i5.5)") it
open(unit=10,file=moviefile,status='unknown')
do iglob = 1,NGLOB
write(10,*) sngl(x(iglob)),sngl(displ(iglob))
@@ -263,13 +267,13 @@
close(10)
endif
- seismogram(itime) = displ(ireceiver)
+ seismogram(it) = displ(ireceiver)
enddo ! end time loop
open(unit=12,file='seismogram',status='unknown')
- do itime = 1,NSTEP
- write(12,*) sngl(dble(itime-1)*DT-hdur),sngl(seismogram(itime))
+ do it = 1,NSTEP
+ write(12,*) sngl(dble(it-1)*DT-hdur),sngl(seismogram(it))
enddo
close(12)
More information about the CIG-COMMITS
mailing list