[cig-commits] [commit] master: Dump counters are no global variables anymore (327e3b0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Aug 25 05:48:30 PDT 2014
Repository : https://github.com/geodynamics/axisem
On branch : master
Link : https://github.com/geodynamics/axisem/compare/d5f53bf9c3506c77131176349d5694c0618b5edc...4fdc3be0b9163702cd293f978f30a53c8313fa17
>---------------------------------------------------------------
commit 327e3b0150c8ec136e044c70e2527f30f36d96a8
Author: Simon Stähler <staehler at geophysik.uni-muenchen.de>
Date: Mon Aug 25 14:01:46 2014 +0200
Dump counters are no global variables anymore
- iseismo, isnap and istrain are now local variables set in sf_time_loop_newmark and symplectic_time_loop and passed down into the dump routines. That should make everything a bit easier to debug.
- Removed several VERY outdated routines used to dump the wavefield into ASCII files. They would probably not have worked anymore anyway.
>---------------------------------------------------------------
327e3b0150c8ec136e044c70e2527f30f36d96a8
SOLVER/data_io.f90 | 9 +--
SOLVER/meshes_io.F90 | 116 ++++++++++++++--------------
SOLVER/nc_routines.F90 | 15 ++--
SOLVER/parameters.F90 | 31 ++------
SOLVER/seismograms.f90 | 70 +++++++++--------
SOLVER/time_evol_wave.F90 | 62 ++++++---------
SOLVER/wavefields_io.f90 | 188 ++++------------------------------------------
7 files changed, 146 insertions(+), 345 deletions(-)
diff --git a/SOLVER/data_io.f90 b/SOLVER/data_io.f90
index 325c1d0..6eea9b3 100644
--- a/SOLVER/data_io.f90
+++ b/SOLVER/data_io.f90
@@ -34,9 +34,6 @@ module data_io
logical :: dump_snaps_glob
logical :: dump_vtk
logical :: dump_xdmf
- logical :: dump_snaps_solflu
- !> N.B. This is not wavefield snapshots, but kernel wavefields. Belongs to nstrain and
- !! istrain
logical :: dump_wavefields
logical :: checkpointing
logical :: diagfiles !< Write diagnostic files (seismograms at antipodes,
@@ -44,9 +41,9 @@ module data_io
logical :: need_fluid_displ
real(kind=dp) :: strain_samp
- integer :: iseismo !< current seismogram sample
- integer :: istrain !< current kernel wavefield sample
- integer :: isnap !< current wavefield sample for movies
+ !integer :: iseismo !< current seismogram sample
+ !integer :: istrain !< current kernel wavefield sample
+ !integer :: isnap !< current wavefield sample for movies
integer :: nseismo !< Number of seismogram samples
integer :: nstrain !< Number of wavefield dumps for kernels
integer :: nsnap !< Number of wavefield snapshots for movies
diff --git a/SOLVER/meshes_io.F90 b/SOLVER/meshes_io.F90
index 9b0f2fc..e80c1bb 100644
--- a/SOLVER/meshes_io.F90
+++ b/SOLVER/meshes_io.F90
@@ -40,8 +40,8 @@ module meshes_io
public :: dump_wavefields_mesh_1d
public :: dump_glob_grid_midpoint
public :: dump_xdmf_grid
- public :: dump_solid_grid
- public :: dump_fluid_grid
+ !public :: dump_solid_grid
+ !public :: dump_fluid_grid
public :: prepare_mesh_memoryvar_vtk
public :: build_kwf_grid
public :: dump_kwf_midpoint_xdmf
@@ -1019,25 +1019,25 @@ end subroutine
!> Dumps the mesh (s,z) [m] in ASCII format as needed to visualize snapshots
!! in the solid region only.
!! Convention for order in the file: First the fluid, then the solid domain.
-subroutine dump_solid_grid(ibeg,iend,jbeg,jend)
-
-
- integer, intent(in) :: ibeg,iend,jbeg,jend
- integer :: iel, ipol,jpol
-
- open(unit=2500+mynum,file=datapath(1:lfdata)//'/solid_grid_'&
- //appmynum//'.dat')
- do iel=1,nel_solid
- do jpol=jbeg,jend
- do ipol=ibeg,iend
- write(2500+mynum,*)scoord(ipol,jpol,ielsolid(iel)), &
- zcoord(ipol,jpol,ielsolid(iel))
- enddo
- enddo
- enddo
- close(2500+mynum)
-
-end subroutine dump_solid_grid
+!subroutine dump_solid_grid(ibeg,iend,jbeg,jend)
+!
+!
+! integer, intent(in) :: ibeg,iend,jbeg,jend
+! integer :: iel, ipol,jpol
+!
+! open(unit=2500+mynum,file=datapath(1:lfdata)//'/solid_grid_'&
+! //appmynum//'.dat')
+! do iel=1,nel_solid
+! do jpol=jbeg,jend
+! do ipol=ibeg,iend
+! write(2500+mynum,*)scoord(ipol,jpol,ielsolid(iel)), &
+! zcoord(ipol,jpol,ielsolid(iel))
+! enddo
+! enddo
+! enddo
+! close(2500+mynum)
+!
+!end subroutine dump_solid_grid
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
@@ -1047,43 +1047,43 @@ end subroutine dump_solid_grid
!! When reading the fluid wavefield, one therefore needs to multiply all
!! components with inv_rho_fluid and the phi component with one/scoord!
!! Convention for order in the file: First the fluid, then the solid domain.
-subroutine dump_fluid_grid(ibeg,iend,jbeg,jend)
-
- use data_pointwise, only : inv_rho_fluid
-
-
- integer, intent(in) :: ibeg,iend,jbeg,jend
- integer :: iel, ipol,jpol
-
- ! When reading the fluid wavefield, one needs to multiply all components
- ! with inv_rho_fluid and the phi component with one/scoord!!
-
- open(unit=2500+mynum,file=datapath(1:lfdata)//&
- '/fluid_grid_'//appmynum//'.dat')
- open(unit=2600+mynum,file=datapath(1:lfdata)//&
- '/inv_rho_scoord_fluid_flusnaps_'&
- //appmynum//'.dat', STATUS="REPLACE")
- do iel=1,nel_fluid
- do jpol=jbeg,jend
- do ipol=ibeg,iend
- write(2500+mynum,*)scoord(ipol,jpol,ielfluid(iel)), &
- zcoord(ipol,jpol,ielfluid(iel))
- if ( axis_fluid(iel) .and. ipol==0 ) then
- ! Axis s=0! write 1 instead of 1/s and then multiply
- ! with the correct factor dsdchi, obtained by L'Hospital's rule
- ! (see routine fluid_snapshot below).
- write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel),one
- else
- write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel), &
- one/scoord(ipol,jpol,ielfluid(iel))
- endif
- enddo
- enddo
- enddo
- close(2500+mynum)
- close(2600+mynum)
-
-end subroutine dump_fluid_grid
+!subroutine dump_fluid_grid(ibeg,iend,jbeg,jend)
+!
+! use data_pointwise, only : inv_rho_fluid
+!
+!
+! integer, intent(in) :: ibeg,iend,jbeg,jend
+! integer :: iel, ipol,jpol
+!
+! ! When reading the fluid wavefield, one needs to multiply all components
+! ! with inv_rho_fluid and the phi component with one/scoord!!
+!
+! open(unit=2500+mynum,file=datapath(1:lfdata)//&
+! '/fluid_grid_'//appmynum//'.dat')
+! open(unit=2600+mynum,file=datapath(1:lfdata)//&
+! '/inv_rho_scoord_fluid_flusnaps_'&
+! //appmynum//'.dat', STATUS="REPLACE")
+! do iel=1,nel_fluid
+! do jpol=jbeg,jend
+! do ipol=ibeg,iend
+! write(2500+mynum,*)scoord(ipol,jpol,ielfluid(iel)), &
+! zcoord(ipol,jpol,ielfluid(iel))
+! if ( axis_fluid(iel) .and. ipol==0 ) then
+! ! Axis s=0! write 1 instead of 1/s and then multiply
+! ! with the correct factor dsdchi, obtained by L'Hospital's rule
+! ! (see routine fluid_snapshot below).
+! write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel),one
+! else
+! write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel), &
+! one/scoord(ipol,jpol,ielfluid(iel))
+! endif
+! enddo
+! enddo
+! enddo
+! close(2500+mynum)
+! close(2600+mynum)
+!
+!end subroutine dump_fluid_grid
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
diff --git a/SOLVER/nc_routines.F90 b/SOLVER/nc_routines.F90
index 4daddc2..16343eb 100644
--- a/SOLVER/nc_routines.F90
+++ b/SOLVER/nc_routines.F90
@@ -530,10 +530,10 @@ end subroutine nc_dump_stf
!-----------------------------------------------------------------------------------------
!> Dump receiver specific stuff, especially displacement and velocity
!! N.B.: Works with global indices.
-subroutine nc_dump_rec(recfield)
+subroutine nc_dump_rec(recfield, iseismo)
use data_mesh, only: num_rec
- use data_io, only: iseismo
real(sp), intent(in), dimension(3,num_rec) :: recfield
+ integer, intent(in) :: iseismo
#ifdef unc
recdumpvar(iseismo,:,:) = 0.0
@@ -613,12 +613,10 @@ end subroutine nc_rec_checkpoint
!----------------------------------------------------------------------------------------
!> Dump stuff along surface
-subroutine nc_dump_surface(surffield, disporvelo)!, nrec, dim2)
- use data_mesh, only: maxind
+subroutine nc_dump_surface(surffield, disporvelo)
- !integer, intent(in) :: nrec, dim2
real(kind=realkind), intent(in), dimension(:,:) :: surffield
- character(len=4), intent(in) :: disporvelo
+ character(len=4), intent(in) :: disporvelo
#ifdef unc
select case(disporvelo)
@@ -1773,15 +1771,16 @@ end subroutine nc_make_snapfile
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
-subroutine nc_dump_snapshot(u, straintrace, curlinplane)
+subroutine nc_dump_snapshot(u, straintrace, curlinplane, isnap)
use data_mesh, only: npoint_plot
- use data_io, only: isnap, nsnap
+ use data_io, only: nsnap
use data_source, only: src_type
real(kind=realkind), dimension(3,npoint_plot), intent(in) :: u
real(kind=realkind), dimension(1,npoint_plot), intent(in) :: straintrace
real(kind=realkind), dimension(1,npoint_plot), intent(in) :: curlinplane
+ integer, intent(in) :: isnap
#ifdef unc
if (src_type(1) == 'monopole') then
diff --git a/SOLVER/parameters.F90 b/SOLVER/parameters.F90
index 3b3b7cd..aa448fd 100644
--- a/SOLVER/parameters.F90
+++ b/SOLVER/parameters.F90
@@ -86,7 +86,6 @@ subroutine readin_parameters
! now pre-set. Most of these are to be considered in the post processing stage now.
sum_seis = .false.
sum_fields = .false.
- dump_snaps_solflu = .false.
num_simul = 1
! netcdf format
@@ -99,7 +98,7 @@ subroutine readin_parameters
write(6,21) datapath, infopath, num_simul, seislength_t, enforced_dt, &
enforced_period, trim(simtype), rec_file_type, &
sum_seis, sum_fields, time_scheme, seis_dt, &
- dump_energy, dump_vtk, dump_snaps_solflu, dump_wavefields, &
+ dump_energy, dump_vtk, dump_wavefields, &
dump_type, ibeg, iend, strain_samp, src_dump_type, make_homo, &
add_hetero, do_mesh_tests, output_format
@@ -180,7 +179,6 @@ subroutine readin_parameters
12x,'Seismogram sampling rate [s]: ',f7.3,/ &
12x,'Dump kin./pot. energy? ',l2,/ &
12x,'Dump global snaps? ',l2,/ &
- 12x,'Dump solid/fluid snaps? ',l2,/ &
12x,'Dump strain? ',l2,/ &
12x,'Wavefield dumping type: ',a12,/ &
12x,'First GLL to save in strains: ',i2,/ &
@@ -202,7 +200,7 @@ subroutine readin_parameters
! Need to decide here since this boolean is needed in def_precomp_terms
need_fluid_displ = .false.
- if (dump_vtk .or. dump_xdmf .or. dump_snaps_solflu .or. dump_energy .or. &
+ if (dump_vtk .or. dump_xdmf .or. dump_energy .or. &
dump_wavefields .and. (dump_type=='fullfields' .or. dump_type=='displ_only')) then
! Need to add this for each new type of wavefield dumping method that
! requires the fluid displacement/velocities
@@ -748,17 +746,6 @@ subroutine check_basic_parameters
14 format(' WARNING: Overriding',a19,' with:',f8.3,' seconds')
- !@ TODO: should this be an ERROR? Do we need dump_snaps_solflu at all given
- ! the newer options?
- if (dump_vtk .and. dump_snaps_solflu) then
- if (lpr) then
- write(6,*)''
- write(6,*)" NOT dumping the same snapshots twice (global AND solid/fluid)"
- write(6,*)'...hence reverting to dumping global snaps only. Sorry.'
- end if
- dump_snaps_solflu = .false.
- endif
-
7 format(04x,a62)
if (verbose > 1) then
@@ -914,7 +901,7 @@ subroutine compute_numerical_parameters
! snapshot output, convert from interval given in seconds to
! incremental time steps
- if (dump_vtk .or. dump_xdmf .or. dump_snaps_solflu .or. dump_memory_vars) then
+ if (dump_vtk .or. dump_xdmf .or. dump_memory_vars) then
snap_it = floor(snap_dt / deltat)
open(unit=2900+mynum, file=datapath(1:lfdata)//'/snap_info.dat'//appmynum)
nsnap = floor(real(niter) / real(snap_it))
@@ -1087,11 +1074,6 @@ subroutine compute_numerical_parameters
endif
endif
- ! Initialize counters for I/O
- istrain = 0
- isnap = 0
- iseismo = 0
-
s_max = zero
! Set some parameters for faster access in time loop
@@ -1350,9 +1332,8 @@ subroutine write_parameters
write(6,12)' Output info path :',trim(infopath)
write(6,19)' Sum wavefields:', sum_fields
write(6,19)' Dump energy :',dump_energy
- write(6,18)' Glob/solflu snaps :',dump_vtk,dump_snaps_solflu
write(6,18)' XDMF VTK :', dump_xdmf
- if (dump_vtk .or. dump_xdmf .or. dump_snaps_solflu) then
+ if (dump_vtk .or. dump_xdmf) then
write(6,11)' snap interval [s] :',snap_dt
write(6,10)' # snaps :',snap_it
endif
@@ -1399,7 +1380,7 @@ subroutine write_parameters
write(55,22)0,'number of strain dumps'
write(55,21)0.,'strain dump sampling rate [s]'
endif
- if (dump_vtk .or. dump_xdmf .or. dump_snaps_solflu) then
+ if (dump_vtk .or. dump_xdmf) then
write(55,22) nsnap,'number of snapshot dumps'
write(55,21)deltat*real(snap_it),'snapshot dump sampling rate [s]'
else
@@ -1473,7 +1454,7 @@ subroutine write_parameters
call nc_write_att_int(0, 'number of strain dumps')
call nc_write_att_dble(0.d0, 'strain dump sampling rate in sec' )
endif
- if (dump_vtk .or. dump_snaps_solflu) then
+ if (dump_vtk) then
call nc_write_att_int(nsnap, 'number of snapshot dumps')
call nc_write_att_real(real(deltat)*real(snap_it), 'snapshot dump sampling rate in sec')
else
diff --git a/SOLVER/seismograms.f90 b/SOLVER/seismograms.f90
index 67233d3..382ae5e 100644
--- a/SOLVER/seismograms.f90
+++ b/SOLVER/seismograms.f90
@@ -895,14 +895,15 @@ end subroutine compute_recfile_seis_bare
!-----------------------------------------------------------------------------
!! Calculate displacement at receiver locations and pass to nc_dump_rec
-subroutine nc_compute_recfile_seis_bare(disp)
+subroutine nc_compute_recfile_seis_bare(disp, iseismo)
use data_source, only : src_type
use nc_routines, only : nc_dump_rec
use data_mesh, only : recfile_el, num_rec, jsurfel
- !use data_mesh
implicit none
real(kind=realkind), intent(in) :: disp(0:,0:,:,:)
+ integer, intent(in) :: iseismo
+
real(kind=realkind) :: disp_rec(3,num_rec)
integer :: i
@@ -929,51 +930,48 @@ subroutine nc_compute_recfile_seis_bare(disp)
enddo
end if !src_type(1)
-call nc_dump_rec(disp_rec)
-
-!deallocate(disp_surf)
+call nc_dump_rec(disp_rec, iseismo)
end subroutine nc_compute_recfile_seis_bare
!=============================================================================
!-----------------------------------------------------------------------------
-subroutine compute_recfile_cmb(velo,grad_sol)
-
- use data_source, only : src_type
- use data_mesh
-
- real(kind=realkind), intent(in) :: velo(0:,0:,:,:)
- real(kind=realkind) :: grad_sol(0:,0:,:,:)
- integer :: i
-
- if (src_type(1)=='monopole') then
- do i=1,num_cmb
- write(200000+i,*)velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1),&
- velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),3)
-
- write(250000+i,*)grad_sol(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1)
- enddo
-
- else
- do i=1,num_cmb
- write(200000+i,*)velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1),&
- velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),2),&
- velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),3)
-
- write(250000+i,*)grad_sol(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1)
-
- enddo
- endif
-
-end subroutine compute_recfile_cmb
+!subroutine compute_recfile_cmb(velo,grad_sol)
+!
+! use data_source, only : src_type
+! use data_mesh
+!
+! real(kind=realkind), intent(in) :: velo(0:,0:,:,:)
+! real(kind=realkind) :: grad_sol(0:,0:,:,:)
+! integer :: i
+!
+! if (src_type(1)=='monopole') then
+! do i=1,num_cmb
+! write(200000+i,*)velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1),&
+! velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),3)
+!
+! write(250000+i,*)grad_sol(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1)
+! enddo
+!
+! else
+! do i=1,num_cmb
+! write(200000+i,*)velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1),&
+! velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),2),&
+! velo(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),3)
+!
+! write(250000+i,*)grad_sol(cmbfile_el(i,2),cmbfile_el(i,3),cmbfile_el(i,1),1)
+!
+! enddo
+! endif
+!
+!end subroutine compute_recfile_cmb
!=============================================================================
!-----------------------------------------------------------------------------
!> 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)
-subroutine compute_surfelem(disp,velo)
+subroutine compute_surfelem(disp, velo)
- use data_io, only : istrain
use data_source, only : src_type
use nc_routines, only : nc_dump_surface
use data_mesh, only : npol, jsurfel, surfelem, maxind
diff --git a/SOLVER/time_evol_wave.F90 b/SOLVER/time_evol_wave.F90
index 2e5ee39..863c0b4 100644
--- a/SOLVER/time_evol_wave.F90
+++ b/SOLVER/time_evol_wave.F90
@@ -171,12 +171,6 @@ subroutine prepare_waves
if (anel_true .and. dump_memory_vars) &
call prepare_mesh_memoryvar_vtk()
- if (dump_snaps_solflu) then
- if (lpr) write(6,*)' dumping solid & fluid grids for snapshots...'
- call dump_solid_grid(ibeg,iend,ibeg,iend)
- if (have_fluid) call dump_fluid_grid(ibeg,iend,ibeg,iend)
- endif
-
! Various seismogram output preparations...
call prepare_seismograms
call open_hyp_epi_equ_anti
@@ -280,7 +274,10 @@ subroutine sf_time_loop_newmark
real(kind=realkind), dimension(0:npol,0:npol,nel_fluid) :: chi, dchi
real(kind=realkind), dimension(0:npol,0:npol,nel_fluid) :: ddchi0, ddchi1
- integer :: iter
+ integer :: iseismo = 0 !< current seismogram sample
+ integer :: istrain = 0 !< current kernel wavefield sample
+ integer :: isnap = 0 !< current wavefield sample for movies
+ integer :: iter
if (lpr) then
write(6,*)
@@ -456,7 +453,7 @@ subroutine sf_time_loop_newmark
acc0(:,:,:,3) = acc1(:,:,:,3)
iclockdump = tick()
- call dump_stuff(iter, disp, velo, chi, dchi, ddchi0, t)
+ call dump_stuff(iter, iseismo, istrain, isnap, disp, velo, chi, dchi, ddchi0, t)
iclockdump = tick(id=iddump, since=iclockdump)
end do ! time loop
@@ -498,7 +495,10 @@ subroutine symplectic_time_loop
! fluid fields
real(kind=realkind), dimension(0:npol,0:npol,nel_fluid) :: chi, dchi, ddchi
- integer :: iter, i
+ integer :: iseismo = 0 !< current seismogram sample
+ integer :: istrain = 0 !< current kernel wavefield sample
+ integer :: isnap = 0 !< current wavefield sample for movies
+ integer :: iter, i
! symplectic stuff
real(kind=dp), allocatable, dimension(:) :: coefd, coeff, coefv, subdt
@@ -667,7 +667,7 @@ subroutine symplectic_time_loop
! ::::::::::::::::::::::::: END SYMPLECTIC SOLVER ::::::::::::::::::::::::::
iclockdump = tick()
- call dump_stuff(iter,disp,velo,chi,dchi,ddchi,t)
+ call dump_stuff(iter, iseismo, istrain, isnap, disp,velo,chi,dchi,ddchi,t)
iclockdump = tick(id=iddump, since=iclockdump)
end do ! time loop
@@ -1015,7 +1015,8 @@ 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
-subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
+subroutine dump_stuff(iter, iseismo, istrain, isnap, &
+ disp, velo, chi, dchi, ddchi, time)
use data_io
use data_mesh
@@ -1024,6 +1025,7 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
use nc_routines, only: nc_rec_checkpoint, nc_dump_strain
integer, intent(in) :: iter
+ integer, intent(inout) :: iseismo, istrain, isnap
real(kind=dp),intent(in) :: time
real(kind=realkind),intent(in) :: disp(0:, 0:, :, :)
@@ -1041,7 +1043,7 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
iseismo = iseismo + 1
if (use_netcdf) then
- call nc_compute_recfile_seis_bare(disp)
+ call nc_compute_recfile_seis_bare(disp, iseismo)
else
call compute_recfile_seis_bare(disp)
endif
@@ -1071,7 +1073,7 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
write(6,*) 'Writing global snap to file: ', isnap
write(6,*)
endif
- call glob_snapshot_midpoint(disp, chi, ibeg, iend, ibeg, iend)
+ call glob_snapshot_midpoint(disp, chi, ibeg, iend, ibeg, iend, isnap)
endif
endif
@@ -1083,23 +1085,10 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
write(6,*)'Writing global xdmf snap to file:',isnap
write(6,*)
endif
- call glob_snapshot_xdmf(disp, chi, time)
+ call glob_snapshot_xdmf(disp, chi, time, isnap)
endif
endif
- !if (dump_snaps_solflu) then
- ! if (mod(iter,snap_it)==0) then
- ! isnap=isnap+1
- ! if (lpr) then
- ! write(6,*)
- ! write(6,*)'Writing solid/fluid snap to file:',isnap
- ! write(6,*)
- ! endif
- ! if (have_fluid) call fluid_snapshot(chi, ibeg, iend, ibeg, iend)
- ! call solid_snapshot(disp, ibeg, iend, ibeg, iend)
- ! endif
- !endif
-
!^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^
! Velocity field and strain tensor wavefields for 3-D kernels^-^-^-^-^
!^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^^-^-^-^-^-^-^-^-^-^-^-^
@@ -1118,8 +1107,6 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
! Any kind of spatial distribution can be dumped, meaning in the long run
! this should be the more effective choice.
!
- ! Currently, 'fullfields' is the hardcoded choice (parameters.f90:110)
- !
! Possible cases in between these dumpsters are considered but not yet
! implemented (e.g. dumping 1/s, inverse fluid density, but compute derivatives
! on-the-fly to warrant grid flexibility).
@@ -1134,7 +1121,7 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
! It starts from
istrain = istrain + 1
- call compute_surfelem(disp,velo)
+ call compute_surfelem(disp, velo)
select case (trim(dump_type))
@@ -1143,15 +1130,15 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
! Minimal permanent storage, minimal run-time memory, minimal CPU time,
! but extensive post-processing (need to compute strain tensor and
! time derivatives, if needed).
- call dump_disp_global(disp, chi) ! displacement globally
+ call dump_disp_global(disp, chi, istrain) ! displacement globally
case ('displ_velo')
! Only dump the 3-comp displacement and velocity fields in solid
! and potential & its derivative in fluid.
! Minimal permanent storage, minimal run-time memory, minimal CPU time,
! but extensive post-processing (need to compute strain tensor).
- call dump_disp(disp, chi) ! displacement in solid, chi in fluid
- call dump_velo_dchi(velo, dchi) ! velocity in solid, dchi in fluid
+ call dump_disp(disp, chi, istrain) ! displacement in solid, chi in fluid
+ call dump_velo_dchi(velo, dchi, istrain) ! velocity in solid, dchi in fluid
case ('fullfields') ! Hardcoded choice
! Compute strain tensor on-the-fly here and dump the 6 components.
@@ -1159,8 +1146,8 @@ subroutine dump_stuff(iter, disp, velo, chi, dchi, ddchi, time)
! Maximal permanent storage, maximal run-time memory, maximal CPU time,
! but no post-processeing necessary as these are the fields that
! constitute density and elastic kernels.
- call compute_strain(disp, chi) ! strain globally
- call dump_velo_global(velo, dchi) ! velocity globally
+ call compute_strain(disp, chi, istrain) ! strain globally
+ call dump_velo_global(velo, dchi, istrain) ! velocity globally
end select
@@ -1227,7 +1214,7 @@ end subroutine dump_stuff
!! scalar Ekk needs to be loaded. Be aware that diuj in the variable names does
!! NOT stand for partial derivatives, but rather the ij component of the
!! strain.
-subroutine compute_strain(u, chi)
+subroutine compute_strain(u, chi, istrain)
use data_pointwise, only: inv_rho_fluid, prefac_inv_s_rho_fluid
use data_source, only: src_type
@@ -1243,6 +1230,7 @@ subroutine compute_strain(u, chi)
real(kind=realkind), intent(in) :: u(0:,0:,:,:)
real(kind=realkind), intent(in) :: chi(0:,0:,:)
+ integer, intent(in) :: istrain
real(kind=realkind) :: grad_sol(0:npol,0:npol,nel_solid,2)
real(kind=realkind) :: buff_solid(0:npol,0:npol,nel_solid)
@@ -1252,7 +1240,7 @@ subroutine compute_strain(u, chi)
character(len=5) :: appisnap
real(kind=realkind), parameter :: two_rk = real(2, kind=realkind)
- call define_io_appendix(appisnap,istrain)
+ call define_io_appendix(appisnap, istrain)
! SSSSSSS Solid region SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS
diff --git a/SOLVER/wavefields_io.f90 b/SOLVER/wavefields_io.f90
index 9efd8eb..b7755ef 100644
--- a/SOLVER/wavefields_io.f90
+++ b/SOLVER/wavefields_io.f90
@@ -35,113 +35,29 @@ module wavefields_io
private
public :: dump_field_1d
- public :: solid_snapshot
public :: glob_snapshot_xdmf
public :: glob_snapshot_midpoint
public :: dump_velo_global
public :: dump_disp_global
public :: dump_disp
public :: dump_velo_dchi
- public :: fluid_snapshot
contains
!-----------------------------------------------------------------------------------------
-!> Dumps the global displacement snapshots [m] in ASCII format
-!! When reading the fluid wavefield, one needs to multiply all
-!! components with inv_rho_fluid and the phi component with one/scoord
-!! as dumped by the corresponding routine dump_glob_grid!
-!! Convention for order in the file: First the fluid, then the solid domain.
-subroutine glob_snapshot(f_sol, chi, ibeg, iend, jbeg, jend)
-
- use data_source, only : src_type
- use pointwise_derivatives, only : axisym_gradient_fluid, dsdf_fluid_axis
- use data_mesh, only : npol, nel_solid, nel_fluid
-
- integer, intent(in) :: ibeg, iend, jbeg, jend
- real(kind=realkind), intent(in) :: f_sol(0:,0:,:,:)
- real(kind=realkind), intent(in) :: chi(0:,0:,:)
-
- real(kind=realkind) :: usz_fluid(0:npol,0:npol,1:nel_fluid,2)
- character(len=4) :: appisnap
- integer :: iel, iidim
- real(kind=realkind) :: dsdchi, prefac
-
- ! When reading the fluid wavefield, one needs to multiply all components
- ! with inv_rho_fluid and the phi component with one/scoord!!
-
- if (src_type(1) == 'monopole') prefac = zero
- if (src_type(1) == 'dipole') prefac = one
- if (src_type(1) == 'quadpole') prefac = two
-
- call define_io_appendix(appisnap, isnap)
-
- open(unit=2500+mynum, file=datapath(1:lfdata)//'/snap_'&
- //appmynum//'_'//appisnap//'.dat')
-
- if (have_fluid) then
- call axisym_gradient_fluid(chi, usz_fluid)
- do iel=1, nel_fluid
-
- if (axis_fluid(iel)) then
- call dsdf_fluid_axis(chi(:,:,iel), iel, 0, dsdchi)
- write(2500+mynum,*) usz_fluid(0,0,iel,1), &
- prefac * dsdchi * chi(0,0,iel), &
- usz_fluid(0,0,iel,2)
- else
- write(2500+mynum,*) usz_fluid(0,0,iel,1), &
- prefac * chi(0,0,iel), &
- usz_fluid(0,0,iel,2)
- endif
-
- write(2500+mynum,*) usz_fluid(npol,0,iel,1), &
- prefac * chi(npol,0,iel), &
- usz_fluid(npol,0,iel,2)
-
- write(2500+mynum,*) usz_fluid(npol,npol,iel,1), &
- prefac * chi(npol,npol,iel), &
- usz_fluid(npol,npol,iel,2)
-
- if ( axis_fluid(iel)) then
- call dsdf_fluid_axis(chi(:,:,iel), iel, npol, dsdchi)
- write(2500+mynum,*) usz_fluid(0,npol,iel,1), &
- prefac * dsdchi * chi(0,npol,iel), &
- usz_fluid(0,npol,iel,2)
- else
- write(2500+mynum,*) usz_fluid(0,npol,iel,1), &
- prefac * chi(0,npol,iel), &
- usz_fluid(0,npol,iel,2)
- endif
-
- enddo
- endif ! have_fluid
-
- do iel=1, nel_solid
- write(2500+mynum,*) (f_sol(0,0,iel,iidim), iidim=1,3)
- write(2500+mynum,*) (f_sol(npol,0,iel,iidim), iidim=1,3)
- write(2500+mynum,*) (f_sol(npol,npol,iel,iidim), iidim=1,3)
- write(2500+mynum,*) (f_sol(0,npol,iel,iidim), iidim=1,3)
- enddo
-
- close(2500+mynum)
-
-end subroutine glob_snapshot
-!-----------------------------------------------------------------------------------------
-
-!-----------------------------------------------------------------------------------------
!> Dumps the global displacement snapshots [m] in binary format
!! When reading the fluid wavefield, one needs to multiply all
!! components with inv_rho_fluid and the phi component with one/scoord
!! as dumped by the corresponding routine dump_glob_grid!
!! Convention for order in the file: First the fluid, then the solid domain.
!! MvD: loop increment npol/2 -> what if npol odd?
-subroutine glob_snapshot_midpoint(f_sol, chi, ibeg, iend, jbeg, jend)
+subroutine glob_snapshot_midpoint(f_sol, chi, ibeg, iend, jbeg, jend, isnap)
use data_source, only : src_type
use pointwise_derivatives, only : axisym_gradient_fluid, dsdf_fluid_axis
use data_mesh, only : npol, nel_solid, nel_fluid
- integer, intent(in) :: ibeg, iend, jbeg, jend
+ integer, intent(in) :: ibeg, iend, jbeg, jend, isnap
real(kind=realkind), intent(in) :: f_sol(0:,0:,:,:)
real(kind=realkind), intent(in) :: chi(0:,0:,:)
@@ -198,7 +114,7 @@ end subroutine glob_snapshot_midpoint
!-----------------------------------------------------------------------------------------
!> Dumps the global displacement snapshots in binary plus XDMF descriptor
-subroutine glob_snapshot_xdmf(f_sol, chi, t)
+subroutine glob_snapshot_xdmf(f_sol, chi, t, isnap)
use data_source, only: src_type
use data_pointwise, only: inv_rho_fluid, prefac_inv_s_rho_fluid
@@ -209,6 +125,7 @@ subroutine glob_snapshot_xdmf(f_sol, chi, t)
real(kind=realkind), intent(in) :: f_sol(0:,0:,:,:)
real(kind=realkind), intent(in) :: chi(0:,0:,:)
real(dp), intent(in) :: t
+ integer, intent(in) :: isnap
character(len=4) :: appisnap
integer :: n_xdmf_fl, n_xdmf_sol
@@ -271,7 +188,7 @@ subroutine glob_snapshot_xdmf(f_sol, chi, t)
i_arr_xdmf, j_arr_xdmf, curlinplane_mask)
! Write variable u to respective file (binary or netcdf)
if (use_netcdf) then
- call nc_dump_snapshot(u, straintrace_mask, curlinplane_mask)
+ call nc_dump_snapshot(u, straintrace_mask, curlinplane_mask, isnap)
else
write(13100) u(1,:)
if (src_type(1) /= 'monopole') write(13101) u(2,:)
@@ -867,89 +784,6 @@ end function
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
-!> Dumps the displacement snapshots [m] in the solid region in ASCII format
-!! Convention for order in the file: First the fluid, then the solid domain.
-subroutine solid_snapshot(f, ibeg, iend, jbeg, jend)
-
- use data_mesh, only: nel_solid
-
- integer, intent(in) :: ibeg, iend, jbeg, jend
- real(kind=realkind), intent(in) :: f(0:,0:,:,:)
- character(len=4) :: appisnap
- integer :: iel, ipol, jpol, idim
-
- call define_io_appendix(appisnap,isnap)
-
- open(unit=3500+mynum, file=datapath(1:lfdata)//'/snap_solid_'&
- //appmynum//'_'//appisnap//'.dat')
-
- do iel=1, nel_solid
- do jpol=ibeg, iend
- do ipol=jbeg, jend
- write(3500+mynum,*) (f(ipol,jpol,iel,idim),idim=1,3)
- enddo
- enddo
- enddo
- close(3500+mynum)
-
-end subroutine solid_snapshot
-!-----------------------------------------------------------------------------------------
-
-!-----------------------------------------------------------------------------------------
-subroutine fluid_snapshot(chi, ibeg, iend, jbeg, jend)
-
- use data_source, only : src_type
- use pointwise_derivatives, only : axisym_gradient_fluid, dsdf_fluid_axis
- use data_mesh, only : npol, nel_fluid
-
- integer, intent(in) :: ibeg, iend, jbeg, jend
- real(kind=realkind), intent(in) :: chi(0:,0:,:)
-
- real(kind=realkind) :: usz_fluid(0:npol,0:npol,1:nel_fluid,2)
- character(len=4) :: appisnap
- integer :: iel, ipol, jpol
- real(kind=realkind) :: dsdchi, prefac
-
- ! When reading the fluid wavefield, one needs to multiply all components
- ! with inv_rho_fluid and the phi component with one/scoord!!
-
- if (src_type(1) == 'monopole') prefac = zero
- if (src_type(1) == 'dipole') prefac = one
- if (src_type(1) == 'quadpole') prefac = two
-
- call axisym_gradient_fluid(chi, usz_fluid)
-
- call define_io_appendix(appisnap, isnap)
-
- open(unit=4500+mynum, file=datapath(1:lfdata)//'/snap_fluid_'&
- //appmynum//'_'//appisnap//'.dat')
-
- do iel=1, nel_fluid
- do jpol=jbeg, jend
- do ipol=ibeg, iend
- if ( axis_fluid(iel) .and. ipol==0 ) then
- call dsdf_fluid_axis(chi(:,:,iel), iel, jpol, dsdchi)
- write(4500+mynum,*) usz_fluid(ipol,jpol,iel,1), &
- !prefac * dsdchi * chi(ipol,jpol,iel), &
- 0., &
- usz_fluid(ipol,jpol,iel,2)
- ! (n.b. up_fl is zero at the axes for all source types, prefac = 0 for
- ! monopole and chi -> 0 for dipole and quadrupole EQ 73-77 in TNM 2007)
- else
- write(4500+mynum,*) usz_fluid(ipol,jpol,iel,1), &
- prefac * chi(ipol,jpol,iel), &
- usz_fluid(ipol,jpol,iel,2)
- endif
- enddo
- enddo
- enddo
-
- close(4500+mynum)
-
-end subroutine fluid_snapshot
-!-----------------------------------------------------------------------------------------
-
-!-----------------------------------------------------------------------------------------
subroutine dump_field_1d(f, filename, appisnap, n)
use data_source, only : have_src, src_dump_type
@@ -990,12 +824,13 @@ end subroutine dump_field_1d
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
-subroutine dump_disp(u, chi)
+subroutine dump_disp(u, chi, istrain)
use data_source, only : src_type,src_dump_type
real(kind=realkind), intent(in) :: u(0:,0:,:,:)
real(kind=realkind), intent(in) :: chi(0:,0:,:)
+ integer, intent(in) :: istrain
integer :: i
character(len=4) :: appisnap
@@ -1035,12 +870,13 @@ end subroutine dump_disp
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
-subroutine dump_velo_dchi(v, dchi)
+subroutine dump_velo_dchi(v, dchi, istrain)
use data_source, only : src_type, src_dump_type
real(kind=realkind),intent(in) :: v(0:,0:,:,:)
real(kind=realkind),intent(in) :: dchi(0:,0:,:)
+ integer, intent(in) :: istrain
integer :: i
character(len=4) :: appisnap
@@ -1080,7 +916,7 @@ end subroutine dump_velo_dchi
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
-subroutine dump_velo_global(v,dchi)
+subroutine dump_velo_global(v, dchi, istrain)
use data_pointwise, only: inv_rho_fluid, prefac_inv_s_rho_fluid
use data_source, only: src_type, src_dump_type
@@ -1089,6 +925,7 @@ subroutine dump_velo_global(v,dchi)
real(kind=realkind), intent(in) :: v(:,:,:,:)
real(kind=realkind), intent(in) :: dchi(:,:,:)
+ integer, intent(in) :: istrain
real(kind=realkind) :: phicomp(0:npol,0:npol,nel_fluid)
integer :: i
@@ -1162,7 +999,7 @@ end subroutine dump_velo_global
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
-subroutine dump_disp_global(u, chi)
+subroutine dump_disp_global(u, chi, istrain)
use data_pointwise, only: inv_rho_fluid, prefac_inv_s_rho_fluid
use data_source, only: src_type, src_dump_type
@@ -1171,6 +1008,7 @@ subroutine dump_disp_global(u, chi)
real(kind=realkind), intent(in) :: u(:,:,:,:)
real(kind=realkind), intent(in) :: chi(:,:,:)
+ integer, intent(in) :: istrain
real(kind=realkind) :: phicomp(0:npol,0:npol,nel_fluid)
integer :: i
More information about the CIG-COMMITS
mailing list