[cig-commits] [commit] master: starting to add a new dump_type: "strain_only" (5f69576)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Oct 13 09:07:27 PDT 2014
Repository : https://github.com/geodynamics/axisem
On branch : master
Link : https://github.com/geodynamics/axisem/compare/ca4b3002407f17f95bbf7013f7622a43ceb4db35...7a0ed4115e0445ccb0219c894eedf4bf157d8689
>---------------------------------------------------------------
commit 5f69576647926acead858fe5ba2a449768dc741b
Author: martinvandriel <martin at vandriel.de>
Date: Mon Oct 13 14:52:33 2014 +0200
starting to add a new dump_type: "strain_only"
this is pointwise strain on a subset of the gll points, with depth and distance filtering.
Useful for instaseis with spatial downsampling
>---------------------------------------------------------------
5f69576647926acead858fe5ba2a449768dc741b
SOLVER/get_mesh.f90 | 8 +-
SOLVER/inparam_advanced.TEMPLATE | 3 +-
SOLVER/meshes_io.F90 | 255 ++++++++++++++++++++++++---------------
SOLVER/nc_routines.F90 | 44 ++++++-
SOLVER/parameters.F90 | 6 +-
SOLVER/time_evol_wave.F90 | 7 +-
6 files changed, 218 insertions(+), 105 deletions(-)
diff --git a/SOLVER/get_mesh.f90 b/SOLVER/get_mesh.f90
index c6a2d04..32e5cb1 100644
--- a/SOLVER/get_mesh.f90
+++ b/SOLVER/get_mesh.f90
@@ -54,7 +54,7 @@ subroutine read_db
use data_comm
use data_proc
use data_time
- use data_io, only : do_anel
+ use data_io, only : do_anel, ibeg, iend, dump_type
use commun, only : barrier, psum, pmax, pmin
use background_models, only : model_is_ani, model_is_anelastic, get_ext_disc, &
override_ext_q
@@ -88,6 +88,12 @@ subroutine read_db
call read_mesh_advanced(1000+mynum)
+ ! npol is now defined, so we can set iend
+ if (trim(dump_type) == 'displ_only') then
+ ibeg = 0
+ iend = npol
+ endif
+
!!!!!!!!!!!! BACKGROUND MODEL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! General numerical input/output parameters
if (verbose > 1) write(69,*)'reading numerical parameters...'
diff --git a/SOLVER/inparam_advanced.TEMPLATE b/SOLVER/inparam_advanced.TEMPLATE
index bbcec50..15fbb98 100644
--- a/SOLVER/inparam_advanced.TEMPLATE
+++ b/SOLVER/inparam_advanced.TEMPLATE
@@ -125,6 +125,7 @@ KERNEL_WAVEFIELDS false
# different modes of dumping wavefields for kernel:
# - fullfields: strain and velocity at a subset of local GLL points (see below)
# - displ_only: displacement at all global points
+# - strain_only: strain at a subset of local GLL points (see below)
KERNEL_DUMPTYPE fullfields
# Samples per period
@@ -137,7 +138,7 @@ KERNEL_SPP 8
KERNEL_SOURCE mask
# GLL points to save, starting and ending GLL point index
-# (ignored for dumptype displ_only)
+# (overwritten with 0 and npol for dumptype displ_only)
KERNEL_IBEG 1
KERNEL_IEND 3
diff --git a/SOLVER/meshes_io.F90 b/SOLVER/meshes_io.F90
index e80c1bb..f83bd1c 100644
--- a/SOLVER/meshes_io.F90
+++ b/SOLVER/meshes_io.F90
@@ -47,6 +47,7 @@ module meshes_io
public :: dump_kwf_midpoint_xdmf
public :: dump_kwf_fem_xdmf
public :: dump_kwf_sem_xdmf
+ public :: dump_kwf_gll_xdmf
contains
@@ -491,6 +492,7 @@ subroutine build_kwf_grid()
use nc_routines, only : set_npoints
use data_mesh
+ use data_io, only: ibeg, iend
integer :: iel, ipol, jpol, ct, ipt, idest
integer, allocatable :: mapping(:)
@@ -554,8 +556,8 @@ subroutine build_kwf_grid()
do iel=1, nel_solid
if (.not. mask_tp_elem(iel)) cycle
- do jpol=0, npol
- do ipol=0, npol
+ do jpol=ibeg, iend
+ do ipol=ibeg, iend
ipt = (iel-1)*(npol+1)**2 + jpol*(npol+1) + ipol + 1
idest = igloc_solid(ipt) + nglob_fluid
@@ -577,8 +579,8 @@ subroutine build_kwf_grid()
do iel=1, nel_fluid
if (.not. mask_tp_elem(iel + nel_solid)) cycle
- do jpol=0, npol
- do ipol=0, npol
+ do jpol=ibeg, iend
+ do ipol=ibeg, iend
ipt = (iel-1)*(npol+1)**2 + jpol*(npol+1) + ipol + 1
idest = igloc_fluid(ipt)
@@ -601,119 +603,121 @@ subroutine build_kwf_grid()
call set_npoints(npoint_kwf)
if (lpr) then
- write(6,*) 'local point number: ', nelem_kwf * (npol + 1)**2
+ write(6,*) 'local point number: ', nelem_kwf * (iend - ibeg + 1)**2
write(6,*) 'after removing duplicates: ', npoint_kwf
write(6,*) 'compression: ', &
real(npoint_kwf) / real(nelem_kwf * (npol + 1)**2)
endif
- allocate(midpoint_mesh_kwf(1:nelem_kwf))
+ if (trim(dump_type) == 'displ_only') then
+ allocate(midpoint_mesh_kwf(1:nelem_kwf))
- if (lpr) write(6,*) ' .... constructing midpoint grid for kwf output'
+ if (lpr) write(6,*) ' .... constructing midpoint grid for kwf output'
- ct = 1
+ ct = 1
- do iel=1, nel_solid
- if (.not. mask_tp_elem(iel)) cycle
- midpoint_mesh_kwf(ct) = mapping_ijel_ikwf(npol/2,npol/2,iel) - 1
- ct = ct + 1
- enddo
+ do iel=1, nel_solid
+ if (.not. mask_tp_elem(iel)) cycle
+ midpoint_mesh_kwf(ct) = mapping_ijel_ikwf(npol/2,npol/2,iel) - 1
+ ct = ct + 1
+ enddo
- do iel=1, nel_fluid
- if (.not. mask_tp_elem(iel + nel_solid)) cycle
- midpoint_mesh_kwf(ct) = mapping_ijel_ikwf(npol/2,npol/2,iel + nel_solid) - 1
- ct = ct + 1
- enddo
+ do iel=1, nel_fluid
+ if (.not. mask_tp_elem(iel + nel_solid)) cycle
+ midpoint_mesh_kwf(ct) = mapping_ijel_ikwf(npol/2,npol/2,iel + nel_solid) - 1
+ ct = ct + 1
+ enddo
- allocate(eltype_kwf(1:nelem_kwf))
+ allocate(eltype_kwf(1:nelem_kwf))
- eltype_kwf(:) = -2
+ eltype_kwf(:) = -2
- ct = 1
+ ct = 1
- do iel=1, nel_solid
- if (.not. mask_tp_elem(iel)) cycle
- if (eltype(ielsolid(iel)) == 'curved') then
- eltype_kwf(ct) = 0
- elseif (eltype(ielsolid(iel)) == 'linear') then
- eltype_kwf(ct) = 1
- elseif (eltype(ielsolid(iel)) == 'semino') then
- eltype_kwf(ct) = 2
- elseif (eltype(ielsolid(iel)) == 'semiso') then
- eltype_kwf(ct) = 3
- else
- eltype_kwf(ct) = -1
- endif
- ct = ct + 1
- enddo
+ do iel=1, nel_solid
+ if (.not. mask_tp_elem(iel)) cycle
+ if (eltype(ielsolid(iel)) == 'curved') then
+ eltype_kwf(ct) = 0
+ elseif (eltype(ielsolid(iel)) == 'linear') then
+ eltype_kwf(ct) = 1
+ elseif (eltype(ielsolid(iel)) == 'semino') then
+ eltype_kwf(ct) = 2
+ elseif (eltype(ielsolid(iel)) == 'semiso') then
+ eltype_kwf(ct) = 3
+ else
+ eltype_kwf(ct) = -1
+ endif
+ ct = ct + 1
+ enddo
- do iel=1, nel_fluid
- if (.not. mask_tp_elem(iel + nel_solid)) cycle
- if (eltype(ielfluid(iel)) == 'curved') then
- eltype_kwf(ct) = 0
- elseif (eltype(ielfluid(iel)) == 'linear') then
- eltype_kwf(ct) = 1
- elseif (eltype(ielfluid(iel)) == 'semino') then
- eltype_kwf(ct) = 2
- elseif (eltype(ielfluid(iel)) == 'semiso') then
- eltype_kwf(ct) = 3
- else
- eltype_kwf(ct) = -1
- endif
- ct = ct + 1
- enddo
+ do iel=1, nel_fluid
+ if (.not. mask_tp_elem(iel + nel_solid)) cycle
+ if (eltype(ielfluid(iel)) == 'curved') then
+ eltype_kwf(ct) = 0
+ elseif (eltype(ielfluid(iel)) == 'linear') then
+ eltype_kwf(ct) = 1
+ elseif (eltype(ielfluid(iel)) == 'semino') then
+ eltype_kwf(ct) = 2
+ elseif (eltype(ielfluid(iel)) == 'semiso') then
+ eltype_kwf(ct) = 3
+ else
+ eltype_kwf(ct) = -1
+ endif
+ ct = ct + 1
+ enddo
- allocate(axis_kwf(1:nelem_kwf))
+ allocate(axis_kwf(1:nelem_kwf))
- axis_kwf(:) = -1
+ axis_kwf(:) = -1
- ct = 1
+ ct = 1
- do iel=1, nel_solid
- if (.not. mask_tp_elem(iel)) cycle
- if (axis_solid(iel)) then
- axis_kwf(ct) = 1
- else
- axis_kwf(ct) = 0
- endif
- ct = ct + 1
- enddo
+ do iel=1, nel_solid
+ if (.not. mask_tp_elem(iel)) cycle
+ if (axis_solid(iel)) then
+ axis_kwf(ct) = 1
+ else
+ axis_kwf(ct) = 0
+ endif
+ ct = ct + 1
+ enddo
- do iel=1, nel_fluid
- if (.not. mask_tp_elem(iel + nel_solid)) cycle
- if (axis_fluid(iel)) then
- axis_kwf(ct) = 1
- else
- axis_kwf(ct) = 0
- endif
- ct = ct + 1
- enddo
+ do iel=1, nel_fluid
+ if (.not. mask_tp_elem(iel + nel_solid)) cycle
+ if (axis_fluid(iel)) then
+ axis_kwf(ct) = 1
+ else
+ axis_kwf(ct) = 0
+ endif
+ ct = ct + 1
+ enddo
- allocate(fem_mesh_kwf(1:4, 1:nelem_kwf))
+ allocate(fem_mesh_kwf(1:4, 1:nelem_kwf))
- if (lpr) write(6,*) ' .... constructing finite element grid for kwf output'
+ if (lpr) write(6,*) ' .... constructing finite element grid for kwf output'
- ct = 1
+ ct = 1
- do iel=1, nel_solid
- if (.not. mask_tp_elem(iel)) cycle
- fem_mesh_kwf(1,ct) = mapping_ijel_ikwf( 0, 0,iel) - 1
- fem_mesh_kwf(2,ct) = mapping_ijel_ikwf(npol, 0,iel) - 1
- fem_mesh_kwf(3,ct) = mapping_ijel_ikwf(npol,npol,iel) - 1
- fem_mesh_kwf(4,ct) = mapping_ijel_ikwf( 0,npol,iel) - 1
- ct = ct + 1
- enddo
+ do iel=1, nel_solid
+ if (.not. mask_tp_elem(iel)) cycle
+ fem_mesh_kwf(1,ct) = mapping_ijel_ikwf( 0, 0,iel) - 1
+ fem_mesh_kwf(2,ct) = mapping_ijel_ikwf(npol, 0,iel) - 1
+ fem_mesh_kwf(3,ct) = mapping_ijel_ikwf(npol,npol,iel) - 1
+ fem_mesh_kwf(4,ct) = mapping_ijel_ikwf( 0,npol,iel) - 1
+ ct = ct + 1
+ enddo
- do iel=1, nel_fluid
- if (.not. mask_tp_elem(iel + nel_solid)) cycle
- fem_mesh_kwf(1,ct) = mapping_ijel_ikwf( 0, 0,iel + nel_solid) - 1
- fem_mesh_kwf(2,ct) = mapping_ijel_ikwf(npol, 0,iel + nel_solid) - 1
- fem_mesh_kwf(3,ct) = mapping_ijel_ikwf(npol,npol,iel + nel_solid) - 1
- fem_mesh_kwf(4,ct) = mapping_ijel_ikwf( 0,npol,iel + nel_solid) - 1
- ct = ct + 1
- enddo
+ do iel=1, nel_fluid
+ if (.not. mask_tp_elem(iel + nel_solid)) cycle
+ fem_mesh_kwf(1,ct) = mapping_ijel_ikwf( 0, 0,iel + nel_solid) - 1
+ fem_mesh_kwf(2,ct) = mapping_ijel_ikwf(npol, 0,iel + nel_solid) - 1
+ fem_mesh_kwf(3,ct) = mapping_ijel_ikwf(npol,npol,iel + nel_solid) - 1
+ fem_mesh_kwf(4,ct) = mapping_ijel_ikwf( 0,npol,iel + nel_solid) - 1
+ ct = ct + 1
+ enddo
+ endif
- allocate(sem_mesh_kwf(0:npol, 0:npol, 1:nelem_kwf))
+ allocate(sem_mesh_kwf(ibeg:iend, ibeg:iend, 1:nelem_kwf))
if (lpr) write(6,*) ' .... constructing spectral element grid for kwf output'
@@ -721,8 +725,8 @@ subroutine build_kwf_grid()
do iel=1, nel_solid
if (.not. mask_tp_elem(iel)) cycle
- do ipol=0, npol
- do jpol=0, npol
+ do ipol=ibeg, iend
+ do jpol=ibeg, iend
sem_mesh_kwf(ipol,jpol,ct) = mapping_ijel_ikwf(ipol,jpol,iel) - 1
enddo
enddo
@@ -731,8 +735,8 @@ subroutine build_kwf_grid()
do iel=1, nel_fluid
if (.not. mask_tp_elem(iel + nel_solid)) cycle
- do ipol=0, npol
- do jpol=0, npol
+ do ipol=ibeg, iend
+ do jpol=ibeg, iend
sem_mesh_kwf(ipol,jpol,ct) = mapping_ijel_ikwf(ipol,jpol,iel + nel_solid) - 1
enddo
enddo
@@ -980,7 +984,8 @@ subroutine dump_kwf_sem_xdmf(filename, npoints, nelem)
open(newunit=iinput_xdmf, file=trim(filename)//'_sem.xdmf')
write(iinput_xdmf, 733) npoints, npoints, trim(filename_np), npoints, trim(filename_np)
- write(iinput_xdmf, 734) nelem, (npol+1)**2, nelem, npol+1, npol+1, trim(filename_np), "'", "'"
+ write(iinput_xdmf, 734) nelem, (iend-ibeg+1)**2, nelem, iend-ibeg+1, iend-ibeg+1, &
+ trim(filename_np), "'", "'"
close(iinput_xdmf)
@@ -1016,6 +1021,54 @@ end subroutine
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
+subroutine dump_kwf_gll_xdmf(filename, npoints)
+ character(len=*), intent(in) :: filename
+ integer, intent(in) :: npoints
+
+ integer :: iinput_xdmf
+ character(len=512) :: filename_np
+
+
+ ! relative filename for xdmf content
+ filename_np = trim(filename(index(filename, '/', back=.true.)+1:))
+
+ ! XML Data
+ open(newunit=iinput_xdmf, file=trim(filename)//'_gll.xdmf')
+ write(iinput_xdmf, 733) npoints, npoints, trim(filename_np), npoints, trim(filename_np)
+
+ write(iinput_xdmf, 734) npoints, "'", "'"
+
+ close(iinput_xdmf)
+
+733 format(&
+ '<?xml version="1.0" ?>',/&
+ '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>',/&
+ '<Xdmf xmlns:xi="http://www.w3.org/2003/XInclude" Version="2.2">',/&
+ '<Domain>',/,/&
+ '<DataItem Name="points" ItemType="Function" Function="join($0, $1)" Dimensions="', i10, ' 2">',/&
+ ' <DataItem Name="points" DataType="Float" Precision="8" Dimensions="', i10, '" Format="HDF">',/&
+ ' ', A, ':/Mesh/mesh_S',/&
+ ' </DataItem>',/&
+ ' <DataItem Name="points" DataType="Float" Precision="8" Dimensions="', i10, '" Format="HDF">',/&
+ ' ', A, ':/Mesh/mesh_Z',/&
+ ' </DataItem>',/&
+ '</DataItem>',/,/)
+
+734 format(&
+ '<Grid Name="grid" GridType="Uniform">',/&
+ ' <Topology TopologyType="Polyvertex" NumberOfElements="',i10'">',/&
+ ' </Topology>',/&
+ ' <Geometry GeometryType="XY">',/&
+ ' <DataItem Reference="/Xdmf/Domain/DataItem[@Name=', A,'points', A,']" />',/&
+ ' </Geometry>',/&
+ '</Grid>',/,/&
+ '</Domain>',/&
+ '</Xdmf>')
+
+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.
@@ -1128,7 +1181,7 @@ subroutine dump_wavefields_mesh_1d
endif
- if (dump_type == 'displ_only') then
+ if (dump_type == 'displ_only' .or. dump_type == 'strain_only') then
call dump_kwf_grid()
else
! compute solid grid
@@ -1188,7 +1241,11 @@ subroutine dump_wavefields_mesh_1d
select case (dump_type)
case ('displ_only')
if (lpr) then
- write(6,*)' strain dump: only global displacement'
+ write(6,*)' strain dump: only elementwise displacement'
+ endif
+ case ('strain_only')
+ if (lpr) then
+ write(6,*)' strain dump: only pointwise strain '
endif
case ('displ_velo')
if (lpr) then
diff --git a/SOLVER/nc_routines.F90 b/SOLVER/nc_routines.F90
index cf44f7e..d446dcd 100644
--- a/SOLVER/nc_routines.F90
+++ b/SOLVER/nc_routines.F90
@@ -777,7 +777,7 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
integer :: iel, ipol, jpol, ct
!print *, 'Processor', mynum,' has been here'
- if (dump_type == 'displ_only') then
+ if (dump_type == 'displ_only' .or. dump_type == 'strain_only') then
allocate(rho1d(npoints))
allocate(lambda1d(npoints))
allocate(mu1d(npoints))
@@ -1052,6 +1052,48 @@ subroutine nc_define_outputfile(nrec, rec_names, rec_th, rec_th_req, rec_ph, rec
nstrain)
end if
+ case ('strain_only')
+ if (src_type(1) == 'monopole') then
+ nvar = 8
+ else
+ nvar = 12
+ end if
+ allocate(varname(nvar))
+ allocate(varnamelist(nvar))
+ allocate(nc_varnamelist(nvar/2))
+ allocate(nc_field_varid(nvar/2))
+
+ if (src_type(1) == 'monopole') then
+ varnamelist = ['strain_dsus_sol', 'strain_dsuz_sol', 'strain_dpup_sol', &
+ 'straintrace_sol', &
+ 'strain_dsus_flu', 'strain_dsuz_flu', 'strain_dpup_flu', &
+ 'straintrace_flu']
+
+ nc_varnamelist = ['strain_dsus', 'strain_dsuz', 'strain_dpup', &
+ 'straintrace']
+ else
+ varnamelist = ['strain_dsus_sol', 'strain_dsuz_sol', 'strain_dpup_sol', &
+ 'strain_dsup_sol', 'strain_dzup_sol', 'straintrace_sol', &
+ 'strain_dsus_flu', 'strain_dsuz_flu', 'strain_dpup_flu', &
+ 'strain_dsup_flu', 'strain_dzup_flu', 'straintrace_flu']
+
+ nc_varnamelist = ['strain_dsus', 'strain_dsuz', 'strain_dpup', &
+ 'strain_dsup', 'strain_dzup', 'straintrace']
+ end if
+
+ npoints = npoint_kwf
+
+ call comm_elem_number(npoints, npoints_global, npoints_myfirst, npoints_mylast)
+ npoint_kwf_global = npoints_global
+
+ npts_sol = npoint_solid_kwf
+ npts_flu = npoint_fluid_kwf
+
+ call comm_elem_number(npts_sol, npts_sol_global, npts_sol_myfirst, npts_sol_mylast)
+ call comm_elem_number(npts_flu, npts_flu_global, npts_flu_myfirst, npts_flu_mylast)
+
+ if (nstrain <= dumpstepsnap) dumpstepsnap = nstrain
+
case ('displ_velo')
write(6,*) 'ERROR: not yet implemented with netcdf'
stop 2
diff --git a/SOLVER/parameters.F90 b/SOLVER/parameters.F90
index 00491e1..99ade14 100644
--- a/SOLVER/parameters.F90
+++ b/SOLVER/parameters.F90
@@ -215,8 +215,9 @@ 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_energy .or. &
- dump_wavefields .and. (dump_type=='fullfields' .or. dump_type=='displ_only')) then
+ if (dump_vtk .or. dump_xdmf .or. dump_energy .or. dump_wavefields .and. &
+ (dump_type=='fullfields' .or. dump_type=='displ_only' &
+ .or. dump_type=='strain_only')) then
! Need to add this for each new type of wavefield dumping method that
! requires the fluid displacement/velocities
need_fluid_displ = .true.
@@ -469,6 +470,7 @@ subroutine read_inparam_advanced
dump_type = to_lower(dump_type)
if (trim(dump_type) /= 'fullfields' .and. &
trim(dump_type) /= 'displ_only' .and. &
+ trim(dump_type) /= 'strain_only' .and. &
trim(dump_type) /= 'displ_velo') then
write(6,*) dump_type
stop 'ERROR: invalid value for KERNEL_DUMPTYPE!'
diff --git a/SOLVER/time_evol_wave.F90 b/SOLVER/time_evol_wave.F90
index a7e0b4e..fd7df7d 100644
--- a/SOLVER/time_evol_wave.F90
+++ b/SOLVER/time_evol_wave.F90
@@ -76,7 +76,8 @@ subroutine prepare_waves
if (rot_src ) call def_rot_matrix
! build mapping to avoid duplicate points at element boundaries
- if (use_netcdf .and. trim(dump_type) == 'displ_only') &
+ if (use_netcdf .and. (trim(dump_type) == 'displ_only' &
+ .or. trim(dump_type) == 'strain_only')) &
call build_kwf_grid()
! Define velocity/density model (velocities in m/s, density in kg/m^3 ) AND
@@ -186,6 +187,10 @@ subroutine prepare_waves
npoint_kwf_global, nelem_kwf_global)
call dump_kwf_sem_xdmf(datapath(1:lfdata)//'/axisem_output.nc4', &
npoint_kwf_global, nelem_kwf_global)
+
+ else if (dump_wavefields .and. dump_type == "strain_only") then
+ call dump_kwf_gll_xdmf(datapath(1:lfdata)//'/axisem_output.nc4', &
+ npoint_kwf_global)
endif
! Need to reload old seismograms and add results
More information about the CIG-COMMITS
mailing list