[cig-commits] [commit] devel, master: adds difference tool to src/tomography (7e26ef7)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:32:47 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 7e26ef7aae2fa8c41a90357607ba5211a4c321a8
Author: daniel peter <peterda at ethz.ch>
Date: Wed Oct 22 17:41:27 2014 +0200
adds difference tool to src/tomography
>---------------------------------------------------------------
7e26ef7aae2fa8c41a90357607ba5211a4c321a8
src/tomography/difference_sem.f90 | 194 ++++++++++++++++++++++++++++++++++++++
src/tomography/rules.mk | 17 ++++
2 files changed, 211 insertions(+)
diff --git a/src/tomography/difference_sem.f90 b/src/tomography/difference_sem.f90
new file mode 100644
index 0000000..09bd35a
--- /dev/null
+++ b/src/tomography/difference_sem.f90
@@ -0,0 +1,194 @@
+! difference_sem
+!
+! this runs in serial, no need to submit as parallel job.
+! takes the difference between proc***.bin from two different input directories
+!
+! usage: ./program_difference_kernel slice_list filename INPUT_dir_1/ INPUT_dir_2/ OUTPUT_DIFF/
+!
+!------------------------------------------------------------------------------------------------
+
+program difference_sem
+
+ use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NX_BATHY,NY_BATHY,IIN,IOUT,MAX_STRING_LEN
+
+ implicit none
+
+ include 'OUTPUT_FILES/values_from_mesher.h'
+
+ integer,parameter :: MAX_NUM_NODES = 2000
+ integer :: node_list(MAX_NUM_NODES)
+
+ character(len=MAX_STRING_LEN) :: arg(6)
+ character(len=MAX_STRING_LEN) :: file1name,file2name
+ character(len=MAX_STRING_LEN) :: input1dir,input2dir
+ character(len=MAX_STRING_LEN) :: outputdir,kernel_name
+ character(len=MAX_STRING_LEN) :: sline
+ character(len=20) :: reg_name
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: sem_data,sem_data_2
+
+ real(kind=CUSTOM_REAL) :: min_all,max_all,min_rel_all,max_rel_all,min,max
+ integer :: num_node,njunk
+ integer :: i, it, iproc, ier
+
+ ! checks arguments
+ do i = 1, 6
+ call get_command_argument(i,arg(i))
+
+ ! usage info
+ if (i < 6 .and. trim(arg(i)) == '') then
+ print *, ' '
+ print *, ' Usage: difference_sem slice_list kernel_name input1_dir/ input2_dir/ output_dir/ '
+ print *, ' '
+ print *, ' with'
+ print *, ' slice_list - text file with slice numbers'
+ print *, ' kernel_name - takes files with this kernel name'
+ print *, ' e.g. "vsv" for proc***_reg1_vsv.bin'
+ print *, ' input1_dir/ - input directory for first files'
+ print *, ' input2_dir/ - input directory for second files'
+ print *, ' output_dir/ - output directory for (first - second) file values'
+ print *, ' '
+ print *, ' possible kernel_names are: '
+ print *, ' "alpha_kernel", "beta_kernel", .., "vsv", "rho_vp", "kappastore", "mustore", etc.'
+ print *
+ print *, ' that are stored in the local directories input1_dir/ and input2_dir/ '
+ print *, ' as real(kind=CUSTOM_REAL) filename(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) in proc***_reg1_filename.bin'
+ print *, ' '
+ stop 'Reenter command line options'
+ endif
+ enddo
+
+ ! get slices id
+ num_node = 0
+ open(unit = 20, file = trim(arg(1)), status = 'old',iostat=ier)
+ if (ier /= 0) then
+ print*,'Error no file: ',trim(arg(1))
+ stop 'Error opening slices file'
+ endif
+
+ ! reads in slices list
+ do while (1 == 1)
+ 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)
+ print *, 'slice list: '
+ print *, node_list(1:num_node)
+ print *, ' '
+
+ ! prefix for crust/mantle region
+ reg_name = 'reg1_'
+
+ ! gets kernel and directory names from argument call
+ kernel_name = trim(reg_name) // trim(arg(2))
+ input1dir = trim(arg(3))
+ input2dir = trim(arg(4))
+ outputdir = trim(arg(5))
+
+ ! loops over slices
+ min_all = 0.0
+ max_all = 0.0
+ min_rel_all = 0.0
+ max_rel_all = 0.0
+
+ do it = 1, num_node
+
+ if(it==1) write(*,*) 'reading files: '
+
+ iproc = node_list(it)
+
+ write(file1name,'(a,i6.6,a)') trim(input1dir)//'/proc',iproc,'_'//trim(kernel_name)//'.bin'
+ write(file2name,'(a,i6.6,a)') trim(input2dir)//'/proc',iproc,'_'//trim(kernel_name)//'.bin'
+
+ ! reads in file from first directory
+ open(IIN,file=trim(file1name),status='old',form='unformatted',iostat=ier)
+ if (ier /= 0 ) then
+ print *,'Error opening file: ',trim(file1name)
+ stop 'Error opening first data file'
+ endif
+ read(IIN) sem_data
+ close(IIN)
+ print *,trim(file1name)
+ print *,' min/max value: ',minval(sem_data),maxval(sem_data)
+
+ ! reads in file from second directory
+ open(IIN,file=trim(file2name),status='old',form='unformatted',iostat=ier)
+ if (ier /= 0 ) then
+ print *,'Error opening file: ',trim(file2name)
+ stop 'Error opening second data file'
+ endif
+ read(IIN) sem_data_2
+ close(IIN)
+ print *,trim(file2name)
+ print *,' min/max value: ',minval(sem_data_2),maxval(sem_data_2)
+
+ ! takes the difference
+ sem_data = sem_data - sem_data_2
+
+ ! stores difference between kernel files
+ if(it==1) write(*,*) 'writing out difference:'
+
+ write(file1name,'(a,i6.6,a)') trim(outputdir)//'/proc',iproc,'_'//trim(kernel_name)//'_diff.bin'
+ open(IOUT,file=trim(file1name),form='unformatted',iostat=ier)
+ if (ier /= 0 ) then
+ print *,'Error opening file: ',trim(file1name)
+ stop 'Error opening output data file'
+ endif
+ write(IOUT) sem_data
+ close(IOUT)
+
+ min = minval(sem_data)
+ max = maxval(sem_data)
+
+ print *,trim(file1name)
+ print *,' min/max value: ',min,max
+
+ if( min < min_all ) min_all = min
+ if( max > max_all ) max_all = max
+
+ ! stores relative difference (k1 - k2)/ k2 with respect to second input file
+ if(it==1) write(*,*) 'writing out relative difference:'
+
+ where( sem_data_2 /= 0.0)
+ sem_data = sem_data / sem_data_2
+ elsewhere
+ sem_data = 0.0
+ endwhere
+
+ write(file1name,'(a,i6.6,a)') trim(outputdir)//'/proc',iproc,'_'//trim(kernel_name)//'_diff_relative.bin'
+ open(IOUT,file=trim(file1name),form='unformatted',iostat=ier)
+ if (ier /= 0 ) then
+ print *,'Error opening file: ',trim(file1name)
+ stop 'Error opening output data file'
+ endif
+
+ write(IOUT) sem_data
+ close(IOUT)
+
+ min = minval(sem_data)
+ max = maxval(sem_data)
+
+ print *,trim(file1name)
+ print *,' min/max relative value: ',minval(sem_data),maxval(sem_data)
+ print *
+
+ if( min < min_rel_all ) min_rel_all = min
+ if( max > max_rel_all ) max_rel_all = max
+
+ enddo
+ write(*,*)
+ write(*,*) 'statistics:'
+ write(*,*) ' total min / max: ',min_all,max_all
+ write(*,*) ' total relative min / max: ',min_rel_all,max_rel_all
+ write(*,*)
+ write(*,*) 'done writing all difference and relative difference files'
+ write(*,*) 'see output directory: ',trim(outputdir)
+ write(*,*)
+
+end program difference_sem
+
+
diff --git a/src/tomography/rules.mk b/src/tomography/rules.mk
index 2dd608a..262f8a0 100644
--- a/src/tomography/rules.mk
+++ b/src/tomography/rules.mk
@@ -32,6 +32,7 @@ tomography_TARGETS = \
$E/xadd_model_tiso \
$E/xadd_model_tiso_cg \
$E/xadd_model_tiso_iso \
+ $E/xdifference_sem \
$E/xinterpolate_model \
$E/xsmooth_sem \
$E/xsum_kernels \
@@ -49,6 +50,7 @@ tomography_OBJECTS = \
$(xadd_model_tiso_OBJECTS) \
$(xadd_model_tiso_cg_OBJECTS) \
$(xadd_model_tiso_iso_OBJECTS) \
+ $(xdifference_sem_OBJECTS) \
$(xinterpolate_model_OBJECTS) \
$(xsmooth_sem_OBJECTS) \
$(xsum_kernels_OBJECTS) \
@@ -59,6 +61,7 @@ tomography_OBJECTS = \
# These files come from the shared directory
tomography_SHARED_OBJECTS = \
$(xadd_model_SHARED_OBJECTS) \
+ $(xdifference_sem_SHARED_OBJECTS) \
$(xinterpolate_model_SHARED_OBJECTS) \
$(xsmooth_sem_SHARED_OBJECTS) \
$(xsum_kernels_SHARED_OBJECTS) \
@@ -186,6 +189,20 @@ ${E}/xconvert_model_file_adios: $(xconvert_model_file_adios_OBJECTS) $(xconvert_
##
+## xdifference_sem
+##
+xdifference_sem_OBJECTS = \
+ $O/difference_sem.tomo.o \
+ $(EMPTY_MACRO)
+
+xdifference_sem_SHARED_OBJECTS = \
+ $O/shared_par.shared_module.o \
+ $(EMPTY_MACRO)
+
+${E}/xdifference_sem: $(xdifference_sem_OBJECTS) $(xdifference_sem_SHARED_OBJECTS)
+ ${FCCOMPILE_CHECK} -o $@ $+
+
+##
## xinterpolate_model
##
xinterpolate_model_OBJECTS = \
More information about the CIG-COMMITS
mailing list