[cig-commits] [commit] master: tidy up (55ce5dc)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Oct 17 05:29:18 PDT 2014
Repository : https://github.com/geodynamics/axisem
On branch : master
Link : https://github.com/geodynamics/axisem/compare/607f803cf074063627513d235f9ed0837fc1dd44...b6457db24acdde4a4e1c08935ae1b22adf87f5bf
>---------------------------------------------------------------
commit 55ce5dc2bd5ab9b1f39d1f84e01aa7d067cbacc2
Author: martinvandriel <martin at vandriel.de>
Date: Wed Oct 15 15:38:16 2014 +0200
tidy up
>---------------------------------------------------------------
55ce5dc2bd5ab9b1f39d1f84e01aa7d067cbacc2
SOLVER/seismograms.f90 | 105 ----------------------------------------
SOLVER/time_evol_wave.F90 | 119 ++++++++++++++--------------------------------
2 files changed, 37 insertions(+), 187 deletions(-)
diff --git a/SOLVER/seismograms.f90 b/SOLVER/seismograms.f90
index 02d56f6..5974b8b 100644
--- a/SOLVER/seismograms.f90
+++ b/SOLVER/seismograms.f90
@@ -999,110 +999,5 @@ subroutine compute_surfelem(disp, velo)
end subroutine compute_surfelem
!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------------------
-!> Save one displacement and velocity trace for each element on the surface
-!! which are both needed for kernels (du and v0 inside the cross-correlation)
-!!
-!!@TODO
-!! MvD: - computes strain in the whole earth, but dumps it only at the surface
-!! - another example copy&paste code -> should go into a proper function
-!! - cannot give the correct result at the axis, i.e. at the antipode and
-!! the epicenter (where 1/s needs treatment with l'hopital)
-!! - what is j? never initatlized....
-!! - sign errors in components 23 and 12 of the strain
-subroutine compute_surfelem_strain(u)
-
- use data_pointwise, only: inv_s_solid
- use data_source, only: src_type
- use pointwise_derivatives, only: axisym_gradient_solid, axisym_gradient_solid_add
- use nc_routines, only: nc_dump_surface
- use data_mesh, only: npol, nel_solid, surfelem, maxind
- real(kind=realkind), intent(in) :: u(0:,0:,:,:)
-
- real(kind=realkind) :: grad_sol(0:npol,0:npol,nel_solid,2)
- real(kind=realkind) :: dumpvar(maxind, 6)
- real(kind=realkind) :: strain(0:npol,nel_solid,6)
-
- integer :: i, jj, j
-
- print *, "You should seriously consider rewriting compute_surfelem_strain()"
- print *, "in seismograms.f90, I am pretty sure it is buggy"
- print *, "If you do not believe me, just go ahead and comment out the stop in"
- print *, " seismograms.f90:1340"
- stop
-
- strain = 0.
-
- if (src_type(1)=='dipole') then
- call axisym_gradient_solid(u(:,:,:,1)+u(:,:,:,2),grad_sol)
- else
- call axisym_gradient_solid(u(:,:,:,1),grad_sol) ! 1: dsus, 2: dzus
- endif
- strain(:,:,1) = grad_sol(npol/2,:,:,1) ! ds us
-
- call axisym_gradient_solid_add(u(:,:,:,3),grad_sol) ! 1:dsuz+dzus,2:dzuz+dsus
-
- ! calculate entire E31 term: (dsuz+dzus)/2
- strain(:,:,4) = grad_sol(npol/2,:,:,1) * real(.5,kind=realkind) ! ds uz
-
- ! Components involving phi....................................................
- if (src_type(1)=='monopole') then
- strain(:,:,2) = inv_s_solid(npol/2,:,:) * u(npol/2,:,:,1) ! dp up
- strain(:,:,3) = grad_sol(npol/2,:,:,2) - strain(:,:,1) ! dz uz
-
- elseif (src_type(1)=='dipole') then
- strain(:,:,2) = real(2.,kind=realkind) * inv_s_solid(npol/2,:,:) * u(npol/2,:,:,1) ! dp up
- strain(:,:,3) = grad_sol(npol/2,:,:,2) - strain(:,:,1) ! dz uz
-
- call axisym_gradient_solid(u(:,:,:,1)-u(:,:,:,2),grad_sol) !1:dsup,2:dzup
- strain(:,:,5) = ( inv_s_solid(npol/2,:,:) * u(npol/2,:,:,2) &
- + real(.5,kind=realkind) * grad_sol(npol/2,:,:,1) ) ! ds up
-
- strain(:,:,6) = real(.5,kind=realkind) * (inv_s_solid(npol/2,:,:) * u(npol/2,:,:,3) &
- + grad_sol(npol/2,:,:,2) ) ! dz up
-
- elseif (src_type(1)=='quadpole') then
- strain(:,:,2) = inv_s_solid(npol/2,:,:) & ! dp up
- * ( u(npol/2,:,:,1) - real(2.,kind=realkind) * u(npol/2,:,:,2))
- strain(:,:,3) = grad_sol(npol/2,:,:,2) - strain(:,:,1) ! dz uz
-
- call axisym_gradient_solid(u(:,:,:,2), grad_sol) ! 1: dsup, 2: dzup
-
- strain(:,:,5) = real(.5,kind=realkind) * ( inv_s_solid(npol/2,:,:) &
- * (real(2.,kind=realkind)* u(npol/2,:,:,1) - u(npol/2,:,:,2)) &
- + grad_sol(npol/2,:,:,1) ) ! ds up
-
- strain(:,:,6) = ( inv_s_solid(npol/2,:,:) *u(npol/2,:,:,3) &
- + real(.5,kind=realkind) *grad_sol(npol/2,:,:,2) ) ! dz up
-
- endif
-
- if (use_netcdf) then
- do i=1, maxind
- dumpvar(i,:) = real(strain(j,surfelem(i),1:6))
- enddo
- call nc_dump_surface(dumpvar(:,1:6), 'stra')
- do i=1, maxind
- dumpvar(i,1:3) = real(u(npol/2,j,surfelem(i),1:3))
- enddo
- call nc_dump_surface(dumpvar(:,1:3), 'srcd')
- end if
-
- if (.not. use_netcdf) then
- do i=1, maxind
- do j=0, npol
- write(60000000+i,20) (real(strain(j,surfelem(i),jj)), jj=1,6)
- write(70000000+i,30) (real(u(npol/2,j,surfelem(i),jj)), jj=1,3)
- enddo
- enddo
-
-20 format(6(1pe11.3))
-30 format(3(1pe11.3))
- end if
-
-end subroutine compute_surfelem_strain
-!-----------------------------------------------------------------------------------------
-
end module seismograms
!=========================================================================================
-
diff --git a/SOLVER/time_evol_wave.F90 b/SOLVER/time_evol_wave.F90
index 7f0b8d7..909c89c 100644
--- a/SOLVER/time_evol_wave.F90
+++ b/SOLVER/time_evol_wave.F90
@@ -19,12 +19,11 @@
! along with AxiSEM. If not, see <http://www.gnu.org/licenses/>.
!
-!========================
+!=========================================================================================
!> Contains all functions for the wave propagation. prepare_waves has to be
!! called beforehand and then time_loop is the only allowed entry point to start
!! wave propagation.
module time_evol_wave
-!========================
use global_parameters
use data_proc
@@ -41,7 +40,7 @@ module time_evol_wave
contains
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Contains all the preliminaries to propagate waves; such as the
!! background model, the stiffness and mass terms, the source and receiver
!! parameters, and preparations for I/O (dumping meshes, opening files).
@@ -215,9 +214,9 @@ subroutine prepare_waves
if (lpr) write(6,*) 'done preparing waves.'
end subroutine prepare_waves
-!=============================================================================
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Entry point into the module time_evol_wave. Calls specific time loop
! !functions, either for newmark or symplectic time scheme.
subroutine time_loop
@@ -235,14 +234,13 @@ subroutine time_loop
iclockold = tick(id=idold, since=iclockold)
end subroutine time_loop
-!=============================================================================
-
+!-----------------------------------------------------------------------------------------
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! T I M E E X T R A P O L A T I O N R O U T I N E S
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> The conventional explicit, acceleration-driven Newmark scheme of 2nd order.
!! (e.g. Chaljub & Valette, 2004). The mass matrix is diagonal; we only store
!! its pre-assembled inverse at the stage of the time loop.
@@ -465,19 +463,19 @@ subroutine sf_time_loop_newmark
end do ! time loop
end subroutine sf_time_loop_newmark
-!=============================================================================
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> SOLVE coupled solid-fluid system of temporal ODEs:
!! M*\dot{u} = -K*u - B*\ddot{\chi} + F (solid)
!! M*\ddot{chi} = -K*\chi - B*u (fluid)
!! using symplectic time integration schemes of 4th, 6th, 8th, 10th order
!!
-!! The time step can be chosen 1.5 times larger than in Newmark, resulting
-!! in CPU times about 2.5 times longer than Newmark, but considerably more
-!! accurate. Consult Ampuero & Nissen-Meyer (2007) for examples of when
-!! this choice might be more appropriate. Generally, for long propagation
-!! distances (say, > 100 wavelengths), it is worthwhile considering this.
+!! The time step can be chosen 1.5 times larger than in Newmark, resulting in CPU
+!! times about 2.5 times longer than Newmark, but considerably more accurate.
+!! Consult Ampuero & Nissen-Meyer (to be submitted 2007) for examples of when this
+!! choice might be more appropriate. Generally, for long propagation distances
+!! (say, > 100 wavelengths), it is worthwhile considering this.
subroutine symplectic_time_loop
use global_parameters
@@ -682,13 +680,13 @@ subroutine symplectic_time_loop
end do ! time loop
end subroutine symplectic_time_loop
-!=============================================================================
+!-----------------------------------------------------------------------------------------
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! E N D O F T I M E E X T R A P O L A T I O N R O U T I N E S
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
subroutine symplectic_coefficients(coefd,coeff,coefv)
use commun, only : barrier,pend
@@ -908,9 +906,9 @@ subroutine symplectic_coefficients(coefd,coeff,coefv)
12 format(' ',i3,' coeffd,coeffv,sub_dt:',3(1pe10.2))
end subroutine symplectic_coefficients
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> coefficients for symmetric compositions of symmetric methods
subroutine SS_scheme(n,a,b,g)
@@ -933,9 +931,9 @@ subroutine SS_scheme(n,a,b,g)
enddo
end subroutine SS_scheme
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Print time step, time, min/max displacement and potential values
!! and stop the simulation if displacements blow up beyond acceptable...
subroutine runtime_info(iter, disp, chi)
@@ -997,9 +995,9 @@ subroutine runtime_info(iter, disp, chi)
endif
end subroutine runtime_info
-!=============================================================================
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Add source term inside source elements only if source time function non-zero
!! and I have the source.
pure subroutine add_source(acc1, stf1)
@@ -1018,9 +1016,9 @@ pure subroutine add_source(acc1, stf1)
endif
end subroutine add_source
-!=============================================================================
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Includes all output action done during the time loop such as
!! various receiver definitions, wavefield snapshots, velocity field & strain
!! tensor for 3-D kernels
@@ -1070,9 +1068,9 @@ subroutine dump_stuff(iter, iseismo, istrain, isnap, &
! Compute kinetic and potential energy globally
if (dump_energy) call energy(disp, velo, dchi, ddchi)
- !^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^
- !^-^-^-^-^-^ Wavefield snapshots-^-^^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^
- !^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^
+ !^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^
+ !^-^-^-^-^-^ Wavefield snapshots-^-^^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^
+ !^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^-^
if (dump_vtk) then
if (mod(iter, snap_it)==0) then
isnap = isnap + 1
@@ -1173,51 +1171,9 @@ subroutine dump_stuff(iter, iseismo, istrain, isnap, &
endif ! dump_wavefields?
end subroutine dump_stuff
-!=============================================================================
+!-----------------------------------------------------------------------------------------
-!!-----------------------------------------------------------------------------
-!!> This is for quick checks and singular computations of the bulk moduli kernels
-!!! hence only the trace of the strain tensor is needed and computed.
-!subroutine dump_velo_straintrace_cmb(u,velo)
-!
-! use data_source, only: src_type
-! use pointwise_derivatives, only: axisym_gradient_solid
-! use pointwise_derivatives, only: f_over_s_solid
-!
-!
-! real(kind=realkind), intent(in) :: u(0:npol,0:npol,nel_solid,3)
-! real(kind=realkind), intent(in) :: velo(0:npol,0:npol,nel_solid,3)
-! real(kind=realkind) :: grad_sol(0:npol,0:npol,nel_solid,3)
-!
-! if (src_type(1)=='dipole') then
-! call axisym_gradient_solid(u(:,:,:,1)+u(:,:,:,2),grad_sol(:,:,:,1:2))
-! else
-! call axisym_gradient_solid(u(:,:,:,1),grad_sol(:,:,:,1:2)) ! dsus, dzus
-! endif
-! call axisym_gradient_solid(u(:,:,:,3),grad_sol(:,:,:,2:3)) ! dsuz, dzuz
-!
-! if (src_type(1)=='monopole') then
-! grad_sol(:,:,:,2)=u(:,:,:,1)
-!
-! elseif (src_type(1)=='dipole') then
-! grad_sol(:,:,:,2)=real(2.,kind=realkind)*u(:,:,:,2)
-!
-! elseif (src_type(1)=='quadpole') then
-! grad_sol(:,:,:,2)=u(:,:,:,1)-real(2.,kind=realkind)*u(:,:,:,2)
-!
-! endif
-!
-! grad_sol(:,:,:,2) = f_over_s_solid(grad_sol(:,:,:,2))
-!
-! ! Slightly dirty, but memory-cheaper: sum of 3 diag strain elements
-! grad_sol(:,:,:,1) = grad_sol(:,:,:,1) + grad_sol(:,:,:,2) + grad_sol(:,:,:,3)
-!
-! call compute_recfile_cmb(velo,grad_sol(:,:,:,1))
-!
-!end subroutine dump_velo_straintrace_cmb
-!!=============================================================================
-
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Compute the full, global strain tensor on-the-fly. Each of 6 (monopole: 4)
!! components is stored separately for solid and fluid domains respectively.
!! The dipole case is transfered to the (s,phi,z) system here.
@@ -1374,10 +1330,9 @@ subroutine compute_strain(u, chi, istrain)
endif ! have_fluid
end subroutine compute_strain
-!=============================================================================
-
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> Computes the kinetic and potential/elastic/stored energy in the solid and
!! fluid subdomains separately. This involves one additional evaluation of
!! the stiffness system (for the velocity vector rather than displacement)
@@ -1488,9 +1443,9 @@ subroutine energy(disp1,vel,dchi1,ddchi)
9 format(4(1pe16.6))
end subroutine energy
-!=============================================================================
+!-----------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> FLUID: solid-fluid boundary term (solid displ.) added to fluid stiffness
!! minus sign in accordance with the definition of bdry_matr
pure subroutine bdry_copy2fluid(uflu,usol)
@@ -1533,9 +1488,9 @@ pure subroutine bdry_copy2fluid(uflu,usol)
endif
end subroutine bdry_copy2fluid
-!============================================================================
+!-----------------------------------------------------------------------------------------
-!----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
!> SOLID: solid-fluid boundary term (fluid potential) added to solid stiffness
!! plus sign in accordance with definition of bdry_matr
pure subroutine bdry_copy2solid(usol,uflu)
@@ -1573,9 +1528,9 @@ pure subroutine bdry_copy2solid(usol,uflu)
end subroutine bdry_copy2solid
-!=============================================================================
+!-----------------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------
-!============================
end module time_evol_wave
-!============================
+!=========================================================================================
More information about the CIG-COMMITS
mailing list