[cig-commits] [commit] devel: updates smoothing, interpolation; makes difference_sem to work in parallel (41a35bb)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Dec 8 08:42:37 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/2cc93c33966a370878432590563633b66f4a1f90...00ae7db05e62851315db5db27dab495a7ec137d3
>---------------------------------------------------------------
commit 41a35bb42931c60e7b14d64bb071eed42236a1ae
Author: daniel peter <peterda at ethz.ch>
Date: Mon Dec 8 15:03:39 2014 +0100
updates smoothing, interpolation; makes difference_sem to work in parallel
>---------------------------------------------------------------
41a35bb42931c60e7b14d64bb071eed42236a1ae
src/tomography/difference_sem.f90 | 252 +++++++++++++++++++++--------------
src/tomography/interpolate_model.F90 | 107 ++++++++++++---
src/tomography/rules.mk | 4 +-
src/tomography/smooth_sem.f90 | 12 +-
4 files changed, 252 insertions(+), 123 deletions(-)
diff --git a/src/tomography/difference_sem.f90 b/src/tomography/difference_sem.f90
index a7aa2c1..70c3367 100644
--- a/src/tomography/difference_sem.f90
+++ b/src/tomography/difference_sem.f90
@@ -43,14 +43,10 @@ program difference_sem
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) :: arg(5)
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
@@ -58,113 +54,163 @@ program difference_sem
real(kind=CUSTOM_REAL) :: min,max,min_rel,max_rel
real(kind=CUSTOM_REAL) :: min_all,max_all,min_rel_all,max_rel_all
- integer :: num_node,njunk
- integer :: i, it, iproc, ier
+ integer :: i,iproc,ier
+ integer :: iregion,ir,irs,ire
+ integer :: nspec
+
+ ! mpi
+ integer :: sizeprocs,myrank
+
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call init_mpi()
+ call world_size(sizeprocs)
+ call world_rank(myrank)
+
+ ! checks compilation setup
+ if (sizeprocs /= NPROCTOT_VAL ) then
+ if (myrank == 0) then
+ print*, 'Error number of processors supposed to run on : ',NPROCTOT_VAL
+ print*, 'Error number of MPI processors actually run on: ',sizeprocs
+ print*
+ print*, 'please rerun with: mpirun -np ',NPROCTOT_VAL,' bin/xdifference_sem .. '
+ endif
+ call exit_MPI(myrank,'Error wrong number of MPI processes')
+ endif
! checks arguments
- do i = 1, 6
+ do i = 1, 5
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'
+ if (i <= 4 .and. trim(arg(i)) == '') then
+ if (myrank == 0) then
+ print *, ' '
+ print *, ' Usage: difference_sem kernel_name input1_dir/ input2_dir/ output_dir/ [region]'
+ print *, ' '
+ print *, ' with'
+ 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 *, ' [region] - optional: if region (1/2/3) is not specified, all 3 regions will be taken,'
+ print *, ' otherwise, only takes region specified'
+ 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 *, ' '
+ endif
+ call synchronize_all()
+ call exit_MPI(myrank,'Reenter command line options')
endif
enddo
- ! get slices id
- num_node = 0
- open(unit = IIN, file = trim(arg(1)), status = 'old',iostat=ier)
- if (ier /= 0) then
- print*,'Error no file: ',trim(arg(1))
- stop 'Error opening slices file'
+ ! get region id
+ if (len_trim(arg(5)) == 0) then
+ iregion = 0
+ else
+ read(arg(5),*) iregion
+ endif
+ if (iregion > 3 .or. iregion < 0) stop 'Error region: region id must be either = 0,1,2,3'
+ if (iregion == 0) then
+ irs = 1
+ ire = 3
+ else
+ irs = iregion
+ ire = irs
endif
-
- ! reads in slices list
- do while (1 == 1)
- read(IIN,'(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(IIN)
- 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))
+ kernel_name = trim(arg(1))
+ input1dir = trim(arg(2))
+ input2dir = trim(arg(3))
+ outputdir = trim(arg(4))
! 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
- ! user output
- write(*,*) 'differencing files: ',it,' out of ',num_node
- iproc = node_list(it)
+ ! user output
+ if (myrank == 0) then
+ write(*,*) 'differencing files: ',sizeprocs,' slices'
+ write(*,*)
+ write(*,*) 'kernel name: ',trim(kernel_name)
+ write(*,*) 'input 1 directory: ',trim(input1dir)
+ write(*,*) 'input 2 directory: ',trim(input2dir)
+ write(*,*) 'output directory : ',trim(outputdir)
+ write(*,*) 'regions: start =', irs, ' to end =', ire
+ write(*,*)
+ endif
- 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'
+ do ir = irs, ire
+ if (myrank == 0) write(*,*) '----------- Region ', ir, '----------------'
+
+ ! initializes
+ sem_data(:,:,:,:) = 0.0_CUSTOM_REAL
+ sem_data_2(:,:,:,:) = 0.0_CUSTOM_REAL
+
+ ! sets number of elements
+ select case (ir)
+ case (1)
+ ! crust/mantle
+ nspec = NSPEC_CRUST_MANTLE
+ case (2)
+ ! outer core
+ nspec = NSPEC_OUTER_CORE
+ case (3)
+ ! inner core
+ nspec = NSPEC_INNER_CORE
+ case default
+ stop 'Error region id not recognized'
+ end select
+
+ ! prefix for (crust/mantle) region
+ write(reg_name,'(a,i1,a)') 'reg',ir,'_'
+
+ ! process file name
+ iproc = myrank
+ write(file1name,'(a,i6.6,a)') trim(input1dir)//'/proc',iproc,'_'//trim(reg_name)//trim(kernel_name)//'.bin'
+ write(file2name,'(a,i6.6,a)') trim(input2dir)//'/proc',iproc,'_'//trim(reg_name)//trim(kernel_name)//'.bin'
! reads in file from first directory
- write(*,*) ' data_1: ',trim(file1name)
+ if(myrank==0) write(*,*) ' data_1: ',trim(file1name)
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
+ read(IIN) sem_data(:,:,:,1:nspec)
close(IIN)
! reads in file from second directory
- write(*,*) ' data_2: ',trim(file2name)
+ if (myrank == 0) write(*,*) ' data_2: ',trim(file2name)
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
+ read(IIN) sem_data_2(:,:,:,1:nspec)
close(IIN)
! user output
- write(*,*) ' min/max data_1 value: ',minval(sem_data),maxval(sem_data)
- write(*,*) ' min/max data_2 value: ',minval(sem_data_2),maxval(sem_data_2)
- write(*,*)
+ if (myrank == 0) then
+ write(*,*) ' min/max data_1 value: ',minval(sem_data(:,:,:,1:nspec)),maxval(sem_data(:,:,:,1:nspec))
+ write(*,*) ' min/max data_2 value: ',minval(sem_data_2(:,:,:,1:nspec)),maxval(sem_data_2(:,:,:,1:nspec))
+ write(*,*)
+ endif
! stores difference between kernel files
- write(*,*) ' difference: (data_1 - data_2)'
+ if (myrank == 0) write(*,*) ' difference: (data_1 - data_2)'
! absolute values
- write(file1name,'(a,i6.6,a)') trim(outputdir)//'/proc',iproc,'_'//trim(kernel_name)//'_diff.bin'
- write(*,*) ' file: ',trim(file1name)
+ write(file1name,'(a,i6.6,a)') trim(outputdir)//'/proc',iproc,'_'//trim(reg_name)//trim(kernel_name)//'_diff.bin'
+ if (myrank == 0) write(*,*) ' file: ',trim(file1name)
open(IOUT,file=trim(file1name),form='unformatted',iostat=ier)
if (ier /= 0 ) then
print *,'Error opening file: ',trim(file1name)
@@ -172,20 +218,20 @@ program difference_sem
endif
! takes the difference
- sem_data = sem_data - sem_data_2
+ sem_data(:,:,:,1:nspec) = sem_data(:,:,:,1:nspec) - sem_data_2(:,:,:,1:nspec)
- write(IOUT) sem_data
+ write(IOUT) sem_data(:,:,:,1:nspec)
close(IOUT)
! min/max
- min = minval(sem_data)
- max = maxval(sem_data)
- if (min < min_all) min_all = min
- if (max > max_all) max_all = max
+ min = minval(sem_data(:,:,:,1:nspec))
+ max = maxval(sem_data(:,:,:,1:nspec))
+ call min_all_cr(min,min_all)
+ call max_all_cr(max,max_all)
! stores relative difference (k1 - k2)/ k2 with respect to second input file
- write(file1name,'(a,i6.6,a)') trim(outputdir)//'/proc',iproc,'_'//trim(kernel_name)//'_diff_relative.bin'
- write(*,*) ' file: ',trim(file1name)
+ write(file1name,'(a,i6.6,a)') trim(outputdir)//'/proc',iproc,'_'//trim(reg_name)//trim(kernel_name)//'_diff_relative.bin'
+ if (myrank == 0) write(*,*) ' file: ',trim(file1name)
open(IOUT,file=trim(file1name),form='unformatted',iostat=ier)
if (ier /= 0 ) then
print *,'Error opening file: ',trim(file1name)
@@ -193,36 +239,44 @@ program difference_sem
endif
! relative difference (k1 - k2)/ k2 with respect to second input file
- where( sem_data_2 /= 0.0)
- sem_data = sem_data / sem_data_2
+ where( sem_data_2(:,:,:,1:nspec) /= 0.0_CUSTOM_REAL)
+ sem_data(:,:,:,1:nspec) = sem_data(:,:,:,1:nspec) / sem_data_2(:,:,:,1:nspec)
elsewhere
- sem_data = 0.0
+ sem_data(:,:,:,1:nspec) = 0.0_CUSTOM_REAL
endwhere
- write(IOUT) sem_data
+ write(IOUT) sem_data(:,:,:,1:nspec)
close(IOUT)
- ! min/max
- min_rel = minval(sem_data)
- max_rel = maxval(sem_data)
- if (min_rel < min_rel_all) min_rel_all = min_rel
- if (max_rel > max_rel_all) max_rel_all = max_rel
+ ! min/max of relative values
+ min_rel = minval(sem_data(:,:,:,1:nspec))
+ max_rel = maxval(sem_data(:,:,:,1:nspec))
+ call min_all_cr(min_rel,min_rel_all)
+ call max_all_cr(max_rel,max_rel_all)
! output
- write(*,*) ' min/max value : ',min,max
- write(*,*) ' min/max relative value: ',min_rel,max_rel
- write(*,*)
- enddo
+ if (myrank == 0) then
+ write(*,*) ' min/max value : ',min,max
+ write(*,*) ' min/max relative value: ',min_rel,max_rel
+ write(*,*)
+ endif
+
+ enddo ! region
! user output
- 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(*,*)
+ if (myrank == 0) then
+ 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(*,*)
+ endif
+
+ ! stop all the MPI processes, and exit
+ call finalize_mpi()
end program difference_sem
diff --git a/src/tomography/interpolate_model.F90 b/src/tomography/interpolate_model.F90
index b5823fe..494067a 100644
--- a/src/tomography/interpolate_model.F90
+++ b/src/tomography/interpolate_model.F90
@@ -105,7 +105,7 @@
logical,parameter :: DO_SEPARATION_TOPO = .true.
! use closest point value in case of large differences
- logical,parameter :: USE_FALLBACK = .false.
+ logical,parameter :: USE_FALLBACK = .true.
!-------------------------------------------------------------------
@@ -1020,8 +1020,7 @@
do iker = 1,nparams
call interpolate(xi,eta,gamma,ispec_selected, &
nspec_max_old,model1(:,:,:,:,iker,rank_selected), &
- val, &
- xigll,yigll,zigll)
+ val,xigll,yigll,zigll)
! sets new model value
model2(i,j,k,ispec,iker) = val
@@ -1350,20 +1349,9 @@
! interpolate model values
do iker = 1,nparams
- call interpolate(xi,eta,gamma,ispec_selected, &
- nspec_max_old,model1(:,:,:,:,iker,rank_selected), &
- val,xigll,yigll,zigll)
-
- ! note: interpolation of values close to the surface or 3D moho encounters problems;
- ! this is a fall-back to the closest point value
- !
- ! uses closest point value if too far off (by more than 5%)
- if (USE_FALLBACK) then
- val_initial = model1(i_selected,j_selected,k_selected,ispec_selected,iker,rank_selected)
- if (abs(val - val_initial ) > abs( 0.05 * val_initial ) ) then
- val = val_initial
- endif
- endif
+ call interpolate_limited(xi,eta,gamma,i_selected,j_selected,k_selected,ispec_selected, &
+ nspec_max_old,model1(:,:,:,:,iker,rank_selected), &
+ val,xigll,yigll,zigll,USE_FALLBACK)
! sets new model value
model2(i,j,k,ispec,iker) = val
@@ -1866,11 +1854,9 @@
!------------------------------------------------------------------------------
!
- subroutine interpolate(xi,eta,gamma, &
- ielem, &
+ subroutine interpolate(xi,eta,gamma,ielem, &
nspec,model, &
- val, &
- xigll,yigll,zigll)
+ val,xigll,yigll,zigll)
use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ
@@ -1912,3 +1898,82 @@
end subroutine interpolate
+!
+!------------------------------------------------------------------------------
+!
+
+ subroutine interpolate_limited(xi,eta,gamma, &
+ i_selected,j_selected,k_selected,ielem, &
+ nspec,model, &
+ val, &
+ xigll,yigll,zigll,USE_FALLBACK)
+
+
+ use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ
+
+ implicit none
+
+ double precision,intent(in):: xi,eta,gamma
+ integer,intent(in):: i_selected,j_selected,k_selected,ielem
+
+ integer,intent(in):: nspec
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: model
+
+ real(kind=CUSTOM_REAL),intent(out) :: val
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX),intent(in) :: xigll
+ double precision, dimension(NGLLY),intent(in) :: yigll
+ double precision, dimension(NGLLZ),intent(in) :: zigll
+
+ logical,intent(in) :: USE_FALLBACK
+
+ ! local parameters
+ double precision :: hxir(NGLLX), hpxir(NGLLX), hetar(NGLLY), hpetar(NGLLY), &
+ hgammar(NGLLZ), hpgammar(NGLLZ)
+ integer:: i,j,k
+ real(kind=CUSTOM_REAL) :: val_initial,val_avg,pert,pert_limit
+
+ ! percentage
+ real(kind=CUSTOM_REAL), parameter :: PERCENTAGE_LIMIT = 0.05
+
+ ! interpolation weights
+ call lagrange_any(xi,NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta,NGLLY,yigll,hetar,hpetar)
+ call lagrange_any(gamma,NGLLZ,zigll,hgammar,hpgammar)
+
+ ! interpolates value
+ val = 0.0_CUSTOM_REAL
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ val = val + hxir(i) * hetar(j) * hgammar(k) * model(i,j,k,ielem)
+ enddo
+ enddo
+ enddo
+
+ ! note: interpolation of values close to the surface or 3D moho encounters problems;
+ ! this is a fall-back to the closest point value
+ !
+ ! uses average/closest point value if too far off (by more than 5%)
+ if (USE_FALLBACK) then
+ ! average value
+ val_avg = sum(model(:,:,:,ielem)) / NGLLX / NGLLY / NGLLZ
+ ! closest point value
+ val_initial = model(i_selected,j_selected,k_selected,ielem)
+
+ ! initial difference
+ pert = abs(val_initial - val_avg)
+
+ ! upper/lower perturbation bound
+ pert_limit = PERCENTAGE_LIMIT * abs(val_avg)
+ if (pert > pert_limit) pert_limit = pert
+
+ ! within a certain percentage range
+ if (abs(val - val_avg ) > pert_limit) then
+ val = val_initial
+ endif
+ endif
+
+ end subroutine interpolate_limited
+
diff --git a/src/tomography/rules.mk b/src/tomography/rules.mk
index d7f40d3..d50feb8 100644
--- a/src/tomography/rules.mk
+++ b/src/tomography/rules.mk
@@ -196,10 +196,12 @@ xdifference_sem_OBJECTS = \
xdifference_sem_SHARED_OBJECTS = \
$O/shared_par.shared_module.o \
+ $O/parallel.sharedmpi.o \
+ $O/exit_mpi.shared.o \
$(EMPTY_MACRO)
${E}/xdifference_sem: $(xdifference_sem_OBJECTS) $(xdifference_sem_SHARED_OBJECTS)
- ${FCCOMPILE_CHECK} -o $@ $+
+ ${MPIFCCOMPILE_CHECK} -o $@ $+
##
## xinterpolate_model
diff --git a/src/tomography/smooth_sem.f90 b/src/tomography/smooth_sem.f90
index 08381b7..8a9f976 100644
--- a/src/tomography/smooth_sem.f90
+++ b/src/tomography/smooth_sem.f90
@@ -226,7 +226,7 @@ program smooth_sem_globe
! average size of a spectral element in km = ...
! e.g. nproc 12x12, nex 192: element_size = 52.122262
if (NCHUNKS_VAL == 6 ) then
- element_size = TWO_PI * dble(4) * R_EARTH_KM / dble(NEX_XI_VAL)
+ element_size = TWO_PI / dble(4) * R_EARTH_KM / dble(NEX_XI_VAL)
else
ANGULAR_WIDTH_XI_RAD = ANGULAR_WIDTH_XI_IN_DEGREES_VAL * DEGREES_TO_RADIANS
ANGULAR_WIDTH_ETA_RAD = ANGULAR_WIDTH_ETA_IN_DEGREES_VAL * DEGREES_TO_RADIANS
@@ -603,6 +603,9 @@ program smooth_sem_globe
endif
endif
+ ! initializes
+ num_elem_local = 0
+
! sets number of elements to loop over
if (.not. DO_BRUTE_FORCE_SEARCH) then
xyz_target(1) = cx0(ispec)
@@ -618,6 +621,11 @@ program smooth_sem_globe
call kdtree_count_nearest_n_neighbors(xyz_target,r_search,num_elem_local)
endif
+ ! checks that at least a single element was choosen
+ if (iproc == myrank) then
+ if (num_elem_local < 1) stop 'Error no local search element found'
+ endif
+
! sets n-search number of nodes
kdtree_search_num_nodes = num_elem_local
@@ -659,8 +667,8 @@ program smooth_sem_globe
deallocate(ispec_flag)
endif
endif
-
else
+ ! brute-force search always loops over whole mesh slice
num_elem_local = NSPEC_AB
endif
More information about the CIG-COMMITS
mailing list