[cig-commits] [commit] devel, master: bug fixes combine_vol_data.F90 to be usable again with adios (0ebb7f4)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:29:07 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 0ebb7f4722762d88e82c79dcff5003cb87c335b1
Author: daniel peter <peterda at ethz.ch>
Date: Thu Aug 21 20:38:02 2014 +0200
bug fixes combine_vol_data.F90 to be usable again with adios
>---------------------------------------------------------------
0ebb7f4722762d88e82c79dcff5003cb87c335b1
Makefile.in | 14 ++-
src/auxiliaries/combine_vol_data.F90 | 160 ++++++++++++++----------
src/auxiliaries/combine_vol_data_adios_impl.f90 | 159 ++++++++++++++++++-----
src/auxiliaries/rules.mk | 16 +--
4 files changed, 242 insertions(+), 107 deletions(-)
diff --git a/Makefile.in b/Makefile.in
index 31ae123..7114287 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -220,6 +220,13 @@ DEFAULT = \
xcreate_movie_GMT_global \
$(EMPTY_MACRO)
+ifeq ($(ADIOS), yes)
+DEFAULT += \
+ xcombine_vol_data_adios \
+ xcombine_vol_data_vtk_adios \
+ $(EMPTY_MACRO)
+endif
+
all: default
default: $(DEFAULT)
@@ -238,18 +245,21 @@ help:
@echo "usage: make [executable]"
@echo ""
@echo "supported executables:"
+ @echo " xcreate_header_file"
@echo " xmeshfem3D"
@echo " xspecfem3D"
- @echo " xcreate_header_file"
@echo " xcombine_vol_data"
@echo " xcombine_vol_data_vtk"
+ifeq ($(ADIOS), yes)
@echo " xcombine_vol_data_adios"
+ @echo " xcombine_vol_data_vtk_adios"
+endif
@echo " xcombine_surf_data"
@echo " xcombine_AVS_DX"
- @echo " xcombine_paraview_strain_data"
@echo " xconvolve_source_timefunction"
@echo " xcreate_movie_AVS_DX"
@echo " xcreate_movie_GMT_global"
+ @echo " xcombine_paraview_strain_data"
@echo ""
.PHONY: all default clean help
diff --git a/src/auxiliaries/combine_vol_data.F90 b/src/auxiliaries/combine_vol_data.F90
index 74494db..848637e 100644
--- a/src/auxiliaries/combine_vol_data.F90
+++ b/src/auxiliaries/combine_vol_data.F90
@@ -25,7 +25,7 @@
!
!=====================================================================
-program combine_vol_data_vtk
+program combine_vol_data
! outputs VTK files (ASCII format)
@@ -36,9 +36,6 @@ program combine_vol_data_vtk
#ifdef ADIOS_INPUT
use adios_read_mod
use combine_vol_data_adios_mod
-!!!!!!!! DK DK removed because this breaks the build system
-!!!!!!!! DK DK use calls to routines in src/shared/parallel.f90 instead
-!!!!!!!! DK DK use mpi
#endif
implicit none
@@ -46,25 +43,25 @@ program combine_vol_data_vtk
include "OUTPUT_FILES/values_from_mesher.h"
integer,parameter :: MAX_NUM_NODES = 2000
- integer :: iregion, ir, irs, ire, ires
- character(len=256) :: sline, arg(7), filename, in_topo_dir, in_file_dir, outdir
- character(len=256) :: prname_topo, prname_file, dimension_file
+ integer :: ir, irs, ire, ires
+ character(len=256) :: arg(7), filename, outdir
character(len=256) :: data_file, topo_file
integer, dimension(MAX_NUM_NODES) :: node_list, nspec, nglob, npoint, nelement
- integer iproc, num_node, i,j,k,ispec, ios, it, di, dj, dk
- integer np, ne, njunk
+ integer :: iproc, num_node, i,j,k,ispec, it, di, dj, dk
+ integer :: np, ne
+ integer :: ier
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: data
real(kind=CUSTOM_REAL),dimension(NGLOB_CRUST_MANTLE) :: xstore, ystore, zstore
- integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE)
+ integer :: ibool(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE)
- integer num_ibool(NGLOB_CRUST_MANTLE)
- logical mask_ibool(NGLOB_CRUST_MANTLE), HIGH_RESOLUTION_MESH
+ integer :: num_ibool(NGLOB_CRUST_MANTLE)
+ logical :: mask_ibool(NGLOB_CRUST_MANTLE), HIGH_RESOLUTION_MESH
- real(kind=CUSTOM_REAL) x, y, z
- real dat
- integer numpoin, iglob, n1, n2, n3, n4, n5, n6, n7, n8
- integer iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
+ real(kind=CUSTOM_REAL) :: x, y, z
+ real :: dat
+ integer :: numpoin, iglob, n1, n2, n3, n4, n5, n6, n7, n8
+ integer :: iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
! instead of taking the first value which appears for a global point, average the values
! if there are more than one GLL points for a global point (points on element corners, edges, faces)
logical,parameter:: AVERAGE_GLOBALPOINTS = .false.
@@ -96,7 +93,6 @@ program combine_vol_data_vtk
#else
!!! .vtk specific !!!!!!!!!!!
character(len=256) :: mesh_file
- integer ier
! global point data
real,dimension(:),allocatable :: total_dat
real,dimension(:,:),allocatable :: total_dat_xyz
@@ -104,32 +100,51 @@ program combine_vol_data_vtk
#endif
#if ADIOS_INPUT
- integer :: sizeprocs, ierr, mpier
+ integer :: sizeprocs,myrank
character(len=256) :: var_name, value_file_name, mesh_file_name
integer(kind=8) :: value_handle, mesh_handle
+#else
+ integer :: iregion,njunk
+ character(len=256) :: prname_topo, prname_file, dimension_file
+ character(len=256) :: in_file_dir, in_topo_dir
+ character(len=256) :: sline
#endif
! starts here---------------------------------------------------------------
+
+ ! dummy initialization to avoid compiler warnings
+ ier = 0
+
+ ! ADIOS mpi initialization
#ifdef ADIOS_INPUT
-!!!!!!!! DK DK removed because this breaks the build system
-!!!!!!!! DK DK use calls to routines in src/shared/parallel.f90 instead
-!!!!!!!! DK DK added this
- stop 'DK DK added this because this part of the code breaks the build system'
-!!!!!!!! DK DK removed because this breaks the build system call MPI_Init(ierr)
-!!!!!!!! DK DK removed because this breaks the build system call MPI_Comm_size(MPI_COMM_WORLD, sizeprocs, ierr)
- print *, sizeprocs, "procs"
+ ! starts mpi
+ call init_mpi()
+ call world_size(sizeprocs)
+ ! checks number of processes
+ ! note: must run as a single process with: mpirun -np 1 ..
if (sizeprocs /= 1) then
- print *, "sequential program. Only mpirun -np 1 ..."
-!!!!!!!! DK DK removed because this breaks the build system call MPI_Abort(MPI_COMM_WORLD, mpier, ierr)
+ call world_rank(myrank)
+ ! usage info
+ if( myrank == 0 ) then
+ print *, "ADIOS requires MPI functionality. However, this program executes as sequential program."
+ print *, "Invalid number of processes used: ", sizeprocs, " procs"
+ print *
+ print *, "Please run: mpirun -np 1 ./bin/xcombine_vol_data_**_adios"
+ endif
+ call abort_mpi()
endif
#endif
+ ! checks array sizes
if (NSPEC_CRUST_MANTLE < NSPEC_OUTER_CORE .or. NSPEC_CRUST_MANTLE < NSPEC_INNER_CORE) &
stop 'This program needs that NSPEC_CRUST_MANTLE > NSPEC_OUTER_CORE and NSPEC_INNER_CORE'
+ ! reads input arguments
#ifndef ADIOS_INPUT
do i = 1, 7
call get_command_argument(i,arg(i))
+
+ ! usage info
if (i < 7 .and. len_trim(arg(i)) == 0) then
print *, ' '
print *, ' Usage: xcombine_vol_data slice_list filename input_topo_dir input_file_dir '
@@ -161,21 +176,23 @@ program combine_vol_data_vtk
! get slices id
num_node = 0
- open(unit = 20, file = trim(arg(1)), status = 'old',iostat = ios)
- if (ios /= 0) then
+ open(unit = 20, file = trim(arg(1)), status = 'old',iostat = ier)
+ if (ier /= 0) then
print*,'no file: ',trim(arg(1))
stop 'Error opening slices file'
endif
do while (1 == 1)
- read(20,'(a)',iostat=ios) sline
- if (ios /= 0) exit
- read(sline,*,iostat=ios) njunk
- if (ios /= 0) exit
+ read(20,'(a)',iostat=ier) sline
+ if (ier /= 0) exit
+ read(sline,*,iostat=ier) njunk
+ if (ier /= 0) exit
num_node = num_node + 1
node_list(num_node) = njunk
enddo
close(20)
+
+ ! output info
print *, 'slice list: '
print *, node_list(1:num_node)
print *, ' '
@@ -191,6 +208,7 @@ program combine_vol_data_vtk
! resolution
read(arg(6),*) ires
#else
+ ! ADIOS input arguments
do i = 1, 7
call get_command_argument(i,arg(i))
enddo
@@ -199,9 +217,11 @@ program combine_vol_data_vtk
outdir, ires, irs, ire)
filename = var_name
#endif
-print *, irs, ire
-!stop
+ ! output info
+ print *, 'regions: start =', irs, ' to end =', ire
+
+ ! sets up loop increments
di = 0
dj = 0
dk = 0
@@ -221,10 +241,12 @@ print *, irs, ire
dj = int((NGLLY-1)/2.0)
dk = int((NGLLZ-1)/2.0)
endif
+
+ ! output info
if( HIGH_RESOLUTION_MESH ) then
- print *, ' high resolution ', HIGH_RESOLUTION_MESH
+ print *, 'using mesh with: high resolution'
else
- print *, ' low resolution ', HIGH_RESOLUTION_MESH
+ print *, 'using mesh with: low resolution'
endif
! sets up ellipticity splines in order to remove ellipticity from point coordinates
@@ -261,9 +283,9 @@ print *, irs, ire
write(prname_file,'(a,i6.6,a,i1,a)') trim(in_file_dir)//'/proc',iproc,'_reg',ir,'_'
dimension_file = trim(prname_topo) //'solver_data.bin'
- open(unit = 27,file = trim(dimension_file),status='old',action='read', iostat = ios, form='unformatted')
- if (ios /= 0) then
- print*,'error ',ios
+ open(unit = 27,file = trim(dimension_file),status='old',action='read', iostat = ier, form='unformatted')
+ if (ier /= 0) then
+ print*,'error ',ier
print*,'file:',trim(dimension_file)
stop 'Error opening file'
endif
@@ -271,8 +293,8 @@ print *, irs, ire
read(27) nglob(it)
close(27)
#else
- call read_scalars_adios_mesh(mesh_handle, iproc, ir, &
- nglob(it), nspec(it))
+ ! adios reading
+ call read_scalars_adios_mesh(mesh_handle, iproc, ir, nglob(it), nspec(it))
#endif
! check
@@ -328,48 +350,56 @@ print *, irs, ire
! write points information
do it = 1, num_node
-
+ ! gets slice id (process id)
iproc = node_list(it)
+ ! initializes data values
data(:,:,:,:) = -1.e10
+ ! output info
print *, ' '
print *, 'Reading slice ', iproc
+
+ ! reads in kernel/data values
#ifndef ADIOS_INPUT
write(prname_topo,'(a,i6.6,a,i1,a)') trim(in_topo_dir)//'/proc',iproc,'_reg',ir,'_'
write(prname_file,'(a,i6.6,a,i1,a)') trim(in_file_dir)//'/proc',iproc,'_reg',ir,'_'
! filename.bin
data_file = trim(prname_file) // trim(filename) // '.bin'
- open(unit = 27,file = trim(data_file),status='old',action='read', iostat = ios,form ='unformatted')
- if (ios /= 0) then
- print*,'error ',ios
+
+ open(unit = 27,file = trim(data_file),status='old',action='read', iostat = ier,form ='unformatted')
+ if (ier /= 0) then
+ print*,'error ',ier
print*,'file:',trim(data_file)
stop 'Error opening file'
endif
- read(27,iostat=ios) data(:,:,:,1:nspec(it))
- if( ios /= 0 ) then
- print*,'read error ',ios
+ read(27,iostat=ier) data(:,:,:,1:nspec(it))
+ if( ier /= 0 ) then
+ print*,'read error ',ier
print*,'file:',trim(data_file)
stop 'error reading data'
endif
close(27)
#else
- call read_values_adios(value_handle, var_name, iproc, ir, nspec(it), data)
+ ! adios reading
+ call read_values_adios(value_handle, var_name, iproc, ir, nspec(it), data)
+ data_file = trim(var_name)
#endif
-
+ ! output info
print *,trim(data_file)
print *,' min/max value: ',minval(data(:,:,:,1:nspec(it))),maxval(data(:,:,:,1:nspec(it)))
print *
! topology file
+ ! reads in mesh coordinates and local-to-global mapping (ibool)
#ifndef ADIOS_INPUT
topo_file = trim(prname_topo) // 'solver_data.bin'
- open(unit = 28,file = trim(topo_file),status='old',action='read', iostat = ios, form='unformatted')
- if (ios /= 0) then
- print*,'error ',ios
- print*,'file:',trim(topo_file)
- stop 'Error opening file'
+ open(unit = 28,file = trim(topo_file),status='old',action='read', iostat = ier, form='unformatted')
+ if (ier /= 0) then
+ print*,'error ',ier
+ print*,'file:',trim(topo_file)
+ stop 'Error opening file'
endif
xstore(:) = 0.0
ystore(:) = 0.0
@@ -384,9 +414,8 @@ print *, irs, ire
if (ir==3) read(28) idoubling_inner_core(1:nspec(it)) ! flag that can indicate fictitious elements
close(28)
#else
- call read_coordinates_adios_mesh(mesh_handle, iproc, ir, &
- nglob(it), nspec(it), &
- xstore, ystore, zstore, ibool)
+ ! adios reading
+ call read_coordinates_adios_mesh(mesh_handle, iproc, ir, nglob(it), nspec(it), xstore, ystore, zstore, ibool)
#endif
print *, trim(topo_file)
@@ -616,8 +645,8 @@ print *, irs, ire
! VTK
! opens unstructured grid file
write(mesh_file,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'.vtk'
- open(IOUT_VTK,file=mesh_file(1:len_trim(mesh_file)),status='unknown',iostat=ios)
- if( ios /= 0 ) stop 'error opening VTK output file'
+ open(IOUT_VTK,file=mesh_file(1:len_trim(mesh_file)),status='unknown',iostat=ier)
+ if( ier /= 0 ) stop 'error opening VTK output file'
write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
write(IOUT_VTK,'(a)') 'material model VTK file'
write(IOUT_VTK,'(a)') 'ASCII'
@@ -678,11 +707,8 @@ print *, irs, ire
#ifdef ADIOS_INPUT
call clean_adios(value_handle, mesh_handle)
-!!!!!!!! DK DK removed because this breaks the build system
-!!!!!!!! DK DK use calls to routines in src/shared/parallel.f90 instead
-!!!!!!!! DK DK added this
- stop 'DK DK added this because this part of the code breaks the build system'
-!!!!!!!! DK DK removed because this breaks the build system call MPI_Finalize(ierr)
+ ! shuts down mpi
+ call finalize_mpi()
#endif
@@ -690,4 +716,4 @@ print *, irs, ire
print *, ' '
-end program combine_vol_data_vtk
+end program combine_vol_data
diff --git a/src/auxiliaries/combine_vol_data_adios_impl.f90 b/src/auxiliaries/combine_vol_data_adios_impl.f90
index 95a1895..a4b88f8 100644
--- a/src/auxiliaries/combine_vol_data_adios_impl.f90
+++ b/src/auxiliaries/combine_vol_data_adios_impl.f90
@@ -26,38 +26,45 @@
!=====================================================================
module combine_vol_data_adios_mod
+
use adios_read_mod
- use mpi
+
implicit none
+
contains
!=============================================================================
!> Print help message.
subroutine print_usage_adios()
+
+ implicit none
print *, 'Usage: '
print *, ' xcombine_data slice_list varname var_file mesh_file ' // &
'output_dir high/low-resolution region'
print *
print *, '* possible varnames are '
- print *, ' rho, rho, kappastore, mustore, alpha_kernel, etc'
+ print *, ' rho, vp, vs, kappastore, mustore, alpha_kl, beta_kl, etc'
print *
print *, ' that are stored in the local directory as ' // &
'real(kind=CUSTOM_REAL) varname(NGLLX,NGLLY,NGLLZ,NSPEC) '
- print *, ' in var_file.bp'
+ print *, ' in a datafile var_file.bp'
print *
- print *, '* mesh_files are used to link variable to the topology'
- print *, '* output_dir indicates where var_name.vtk will be written'
+ print *, '* mesh_files: are used to link variable to the topology (e.g. DATABASES_MPI/solver_data.bp)'
+ print *, '* output_dir: indicates where var_name.vtk will be written'
print *, '* give 0 for low resolution and 1 for high resolution'
print *
stop ' Reenter command line options'
+
end subroutine print_usage_adios
!=============================================================================
!> Interpret command line arguments
+
subroutine read_args_adios(arg, MAX_NUM_NODES, node_list, num_node, &
var_name, value_file_name, mesh_file_name, &
outdir, ires, irs, ire)
+
implicit none
! Arguments
character(len=*), intent(in) :: arg(:)
@@ -70,8 +77,7 @@ subroutine read_args_adios(arg, MAX_NUM_NODES, node_list, num_node, &
character(len=256) :: sline
integer :: ios, njunk, iregion
- if ((command_argument_count() == 6) &
- .or. (command_argument_count() == 7)) then
+ if ((command_argument_count() == 6) .or. (command_argument_count() == 7)) then
num_node = 0
open(unit = 20, file = trim(arg(1)), status = 'unknown',iostat = ios)
if (ios /= 0) then
@@ -84,8 +90,7 @@ subroutine read_args_adios(arg, MAX_NUM_NODES, node_list, num_node, &
read(sline,*,iostat=ios) njunk
if (ios /= 0) exit
num_node = num_node + 1
- if( num_node > MAX_NUM_NODES ) &
- stop 'error number of slices exceeds MAX_NUM_NODES...'
+ if( num_node > MAX_NUM_NODES ) stop 'error number of slices exceeds MAX_NUM_NODES...'
node_list(num_node) = njunk
enddo
close(20)
@@ -116,27 +121,36 @@ end subroutine read_args_adios
!=============================================================================
!> Open ADIOS value and mesh files, read mode
-subroutine init_adios(value_file_name, mesh_file_name, &
- value_handle, mesh_handle)
+
+subroutine init_adios(value_file_name, mesh_file_name, value_handle, mesh_handle)
+
implicit none
! Parameters
character(len=*), intent(in) :: value_file_name, mesh_file_name
integer(kind=8), intent(out) :: value_handle, mesh_handle
! Variables
integer :: ier
+ integer :: comm
- call adios_read_init_method(ADIOS_READ_METHOD_BP, MPI_COMM_WORLD, "verbose=1", ier)
+ call world_get_comm(comm)
+
+ call adios_read_init_method(ADIOS_READ_METHOD_BP, comm, "verbose=1", ier)
+ if( ier /= 0 ) stop 'error in adios_read_init_method()'
+
+ call adios_read_open_file(mesh_handle, trim(mesh_file_name), 0, comm, ier)
+ if( ier /= 0 ) stop 'error opening adios mesh file with adios_read_open_file()'
+
+ call adios_read_open_file(value_handle, trim(value_file_name), 0, comm, ier)
+ if( ier /= 0 ) stop 'error opening value file with adios_read_open_file()'
- call adios_read_open_file(mesh_handle, trim(mesh_file_name), 0, &
- MPI_COMM_WORLD, ier)
- call adios_read_open_file(value_handle, trim(value_file_name), 0, &
- MPI_COMM_WORLD, ier)
end subroutine init_adios
!=============================================================================
!> Open ADIOS value and mesh files, read mode
+
subroutine clean_adios(value_handle, mesh_handle)
+
implicit none
! Parameters
integer(kind=8), intent(in) :: value_handle, mesh_handle
@@ -150,7 +164,9 @@ end subroutine clean_adios
!=============================================================================
+
subroutine read_scalars_adios_mesh(mesh_handle, iproc, ir, nglob, nspec)
+
implicit none
! Parameters
integer(kind=8), intent(in) :: mesh_handle
@@ -164,15 +180,26 @@ subroutine read_scalars_adios_mesh(mesh_handle, iproc, ir, nglob, nspec)
write(reg_name, '(a,i1)') trim("reg"), ir
call adios_selection_writeblock(sel, iproc)
- call adios_schedule_read(mesh_handle, sel, trim(reg_name) // "/nglob", &
- 0, 1, nglob, ier)
- call adios_schedule_read(mesh_handle, sel, trim(reg_name) // "/nspec", &
- 0, 1, nspec, ier)
+ call adios_schedule_read(mesh_handle, sel, trim(reg_name) // "/nglob", 0, 1, nglob, ier)
+ if( ier /= 0 ) then
+ print*
+ print* ,'error adios: could not read parameter: ',trim(reg_name) // "/nglob"
+ print* ,' please check if your input mesh file is correct...'
+ print*
+ stop 'error adios adios_schedule_read() for nglob'
+ endif
+
+ call adios_schedule_read(mesh_handle, sel, trim(reg_name) // "/nspec", 0, 1, nspec, ier)
+ if( ier /= 0 ) stop 'error adios adios_schedule_read() for nspec'
+
call adios_perform_reads(mesh_handle, ier)
+ if( ier /= 0 ) stop 'error adios: perform reading nglob and nspec failed'
+
end subroutine read_scalars_adios_mesh
!=============================================================================
+
subroutine read_coordinates_adios_mesh(mesh_handle, iproc, ir, nglob, nspec, &
xstore, ystore, zstore, ibool)
use constants
@@ -195,10 +222,22 @@ subroutine read_coordinates_adios_mesh(mesh_handle, iproc, ir, nglob, nspec, &
call adios_schedule_read(mesh_handle, sel_scalar, &
trim(reg_name) // "ibool/offset", &
0, 1, offset_ibool, ier)
+ if( ier /= 0 ) then
+ print*
+ print* ,'error adios: could not read parameter: ',trim(reg_name) // "ibool/offset"
+ print* ,' please check if your input mesh file is correct...'
+ print*
+ stop 'error adios adios_schedule_read() for ibool/offset'
+ endif
+
+
call adios_schedule_read(mesh_handle, sel_scalar, &
trim(reg_name) // "x_global/offset", &
0, 1, offset_coord, ier)
+ if( ier /= 0 ) stop 'error adios: reading x_global/offset'
+
call adios_perform_reads(mesh_handle, ier)
+ if( ier /= 0 ) stop 'error adios: perform reading mesh file offsets failed'
start(1) = offset_ibool
count_ad(1) = NGLLX * NGLLY * NGLLZ * nspec
@@ -224,8 +263,10 @@ end subroutine read_coordinates_adios_mesh
!=============================================================================
-subroutine read_values_adios(value_handle, var_name, iproc, ir, &
- nspec, data)
+!> reads in data from ADIOS value file
+
+subroutine read_values_adios(value_handle, var_name, iproc, ir, nspec, data)
+
use constants
implicit none
@@ -235,26 +276,84 @@ subroutine read_values_adios(value_handle, var_name, iproc, ir, &
integer, intent(in) :: iproc, ir, nspec
real(kind=CUSTOM_REAL), dimension(:,:,:,:), intent(inout) :: data
! Variables
- character(len=256) :: reg_name
integer(kind=8), dimension(1) :: start, count_ad
integer(kind=8) :: sel
integer :: offset, ier
+ character(len=256) :: data_name
+ character(len=8) :: reg_name
+ logical :: is_kernel
+
+ ! note: we can either visualize
+ ! wavespeed arrays (in DATABASES_MPI/solver_meshfiles.bp)
+ ! (e.g. rho,vp,vs,vph,vpv,vsh,..)
+ ! or
+ ! sensitivity kernels (in OUTPUT_FILES/kernels.bp)
+ ! (e.g. rho_kl,alpha_kl,beta_kl,alphah_kl,alphav_kl,betah_kl,..)
+ !
+ ! unfortunately, they have different naming conventions for different Earth regions:
+ ! - wavespeed arrays: reg1/vp/, reg2/vp/, reg3/vp/, ..
+ ! - kernels: alpha_kl_crust_mantle/, alpha_kl_outer_core/, alpha_kl_inner_core/,..
+ ! here we try to estimate the type by the ending of the variable name given by the user,
+ ! i.e. if the ending is "_kl" we assume it is a kernel name
+ !
+ is_kernel = .false.
+ if( len_trim(var_name) > 3 ) then
+ if( var_name(len_trim(var_name)-2:len_trim(var_name)) == '_kl' ) then
+ is_kernel = .true.
+ endif
+ endif
- write(reg_name, '(a,i1, a)') trim("reg"), ir, "/"
+ ! determines full data array name
+ if( is_kernel ) then
+ ! for kernel name: alpha_kl, .. , **_kl only:
+ ! adds region name to get kernel_name
+ ! for example: var_name = "alpha_kl" -> alpha_kl_crust_mantle
+ !
+ ! note: this must match the naming convention used
+ ! in file save_kernels_adios.F90
+ select case( ir )
+ case( IREGION_CRUST_MANTLE )
+ data_name = trim(var_name) // "_crust_mantle"
+ case( IREGION_OUTER_CORE )
+ data_name = trim(var_name) // "_outer_core"
+ case( IREGION_INNER_CORE )
+ data_name = trim(var_name) // "_inner_core"
+ case default
+ stop 'error wrong region code in read_values_adios() routine'
+ end select
+ else
+ ! for wavespeed name: rho,vp,..
+ ! adds region name: var_name = "rho" -> reg1/rho
+ write(reg_name, '(a,i1, a)') trim("reg"), ir, "/"
+ data_name = trim(reg_name) // trim(var_name)
+ endif
+ ! reads in data offset
call adios_selection_writeblock(sel, iproc)
- call adios_schedule_read(value_handle, sel, &
- trim(reg_name) // trim(var_name) // "/offset", &
- 0, 1, offset, ier)
+
+ call adios_schedule_read(value_handle, sel, trim(data_name) // "/offset", 0, 1, offset, ier)
+ if( ier /= 0 ) then
+ print*
+ print* ,'error adios: could not read parameter: ',trim(data_name) // "/offset"
+ print* ,' please check if your input data file is correct...'
+ print*
+ stop 'error adios adios_schedule_read() for data offset'
+ endif
+
call adios_perform_reads(value_handle, ier)
+ if( ier /= 0 ) stop 'error adios: perform reading data offset failed'
+ ! reads in data array
start(1) = offset
count_ad(1) = NGLLX * NGLLY * NGLLZ * nspec
call adios_selection_boundingbox (sel , 1, start, count_ad)
- call adios_schedule_read(value_handle, sel, &
- trim(reg_name) // trim(var_name) // "/array", 0, 1,&
- data, ier)
+
+ call adios_schedule_read(value_handle, sel, trim(data_name) // "/array", 0, 1, data, ier)
+ if( ier /= 0 ) stop 'error adios reading data array'
+
call adios_perform_reads(value_handle, ier)
+ if( ier /= 0 ) stop 'error adios perform reading of data array'
+
end subroutine read_values_adios
end module combine_vol_data_adios_mod
diff --git a/src/auxiliaries/rules.mk b/src/auxiliaries/rules.mk
index 9f3df67..cbdde73 100644
--- a/src/auxiliaries/rules.mk
+++ b/src/auxiliaries/rules.mk
@@ -104,10 +104,10 @@ ${E}/xcombine_paraview_strain_data: $(auxiliaries_SHARED_OBJECTS) $O/combine_par
${E}/xcombine_vol_data: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data.auxsolver.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
${FCCOMPILE_CHECK} -o $@ $+
-${E}/xcombine_vol_data_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
+${E}/xcombine_vol_data_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o $O/parallel.sharedmpi.o
${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
-${E}/xcombine_vol_data_vtk_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
+${E}/xcombine_vol_data_vtk_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o $O/parallel.sharedmpi.o
${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
${E}/xcombine_vol_data_vtk: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data.auxsolver_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
@@ -151,8 +151,8 @@ $O/%.aux.o: $S/%.f90 $O/shared_par.shared_module.o ${OUTPUT}/values_from_mesher.
$O/%.auxsolver.o: $S/%.f90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o
${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $<
-$O/%.auxmpi.o: $S/%.f90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o
- ${MPIFCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $<
+$O/%.auxmpi.o: $S/%.f90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o $O/parallel.sharedmpi.o
+ ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $<
$O/%.auxsolver.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o
${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $<
@@ -160,8 +160,8 @@ $O/%.auxsolver.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_m
$O/%.auxsolver_vtk.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o
${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< $(FC_DEFINE)USE_VTK_INSTEAD_OF_MESH
-$O/%.auxadios.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o
- ${MPIFCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< $(FC_DEFINE)ADIOS_INPUT
+$O/%.auxadios.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o $O/parallel.sharedmpi.o
+ ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< $(FC_DEFINE)ADIOS_INPUT
-$O/%.auxadios_vtk.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o
- ${MPIFCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< $(FC_DEFINE)ADIOS_INPUT $(FC_DEFINE)USE_VTK_INSTEAD_OF_MESH
+$O/%.auxadios_vtk.o: $S/%.F90 ${OUTPUT}/values_from_mesher.h $O/shared_par.shared_module.o $O/parallel.sharedmpi.o
+ ${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $< $(FC_DEFINE)ADIOS_INPUT $(FC_DEFINE)USE_VTK_INSTEAD_OF_MESH
More information about the CIG-COMMITS
mailing list