[cig-commits] r18639 - in seismo/3D/SPECFEM3D_GLOBE/trunk: setup src/meshfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Jun 15 15:59:58 PDT 2011
Author: dkomati1
Date: 2011-06-15 15:59:57 -0700 (Wed, 15 Jun 2011)
New Revision: 18639
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90
Log:
committed Mila's modifications to merge by GPU mesh coloring routines with the official mesher of SPECFEM3D_GLOBE, in order to maintain a single code and prepare for the future official GPU version of the SPECFEM3D_GLOBE solver
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2011-06-15 22:30:12 UTC (rev 18638)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2011-06-15 22:59:57 UTC (rev 18639)
@@ -51,6 +51,20 @@
integer, parameter :: ELEMENTS_NONBLOCKING_CM_IC = 1500
integer, parameter :: ELEMENTS_NONBLOCKING_OC = 3000
+!*********************************************************************************************************
+! added these parameters for the future GPU version of the solver with mesh coloring
+
+! sort outer elements first and then inner elements
+! in order to use non blocking MPI to overlap communications
+ logical, parameter :: SORT_MESH_INNER_OUTER = .true.
+
+! add mesh coloring for the GPU + MPI implementation
+ logical, parameter :: USE_MESH_COLORING_GPU = .false.
+ integer, parameter :: MAX_NUMBER_OF_COLORS = 10000
+ integer, parameter :: NGNOD_HEXAHEDRA = 8
+
+!*********************************************************************************************************
+
! if files on a local path on each node are also seen as global with same path
! set to .true. typically on a shared-memory machine with a common file system
! set to .false. typically on a cluster of nodes with local disks
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/Makefile.in 2011-06-15 22:30:12 UTC (rev 18638)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/Makefile.in 2011-06-15 22:59:57 UTC (rev 18639)
@@ -105,6 +105,7 @@
$O/create_doubling_elements.o \
$O/create_mass_matrices.o \
$O/create_regions_mesh.o \
+ $O/get_perm_color.o \
$O/create_regular_elements.o \
$O/define_superbrick.o \
$O/euler_angles.o \
@@ -346,8 +347,11 @@
${FCCOMPILE_CHECK} -c -o $O/create_mass_matrices.o ${FCFLAGS_f90} $S/create_mass_matrices.f90
$O/create_regions_mesh.o: ${SETUP}/constants.h $S/create_regions_mesh.f90
- ${FCCOMPILE_CHECK} -c -o $O/create_regions_mesh.o ${FCFLAGS_f90} $S/create_regions_mesh.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/create_regions_mesh.o ${FCFLAGS_f90} $S/create_regions_mesh.f90
+$O/get_perm_color.o: ${SETUP}/constants.h $S/get_perm_color.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_perm_color.o ${FCFLAGS_f90} $S/get_perm_color.f90
+
$O/create_regular_elements.o: ${SETUP}/constants.h $S/create_regular_elements.f90
${FCCOMPILE_CHECK} -c -o $O/create_regular_elements.o ${FCFLAGS_f90} $S/create_regular_elements.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90 2011-06-15 22:30:12 UTC (rev 18638)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90 2011-06-15 22:59:57 UTC (rev 18639)
@@ -40,7 +40,7 @@
R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,&
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
ner,ratio_sampling_array,doubling_index,r_bottom,r_top, &
- this_region_has_a_doubling,ipass,ratio_divide_central_cube,&
+ this_region_has_a_doubling,ipass,ratio_divide_central_cube, &
CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,offset_proc_xi,offset_proc_eta)
! creates the different regions of the mesh
@@ -49,6 +49,15 @@
implicit none
+!****************************************************************************************************
+! Mila
+
+! include "constants.h"
+! standard include of the MPI library
+ include 'mpif.h'
+
+!****************************************************************************************************
+
! this to cut the doubling brick
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_XI_FACE,NSPEC2D_ETA_FACE
@@ -216,6 +225,20 @@
logical :: ACTUALLY_STORE_ARRAYS
+!****************************************************************************************************
+! Mila
+
+! added for color permutation
+ integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer
+ integer, dimension(:), allocatable :: perm
+ integer, dimension(:), allocatable :: first_elem_number_in_this_color
+ integer, dimension(:), allocatable :: num_of_elems_in_this_color
+
+ integer :: icolor,ispec_counter
+ integer :: nspec_outer_min_global,nspec_outer_max_global
+
+!****************************************************************************************************
+
! Boundary Mesh
integer NSPEC2D_MOHO,NSPEC2D_400,NSPEC2D_670,nex_eta_moho
integer, dimension(:), allocatable :: ibelm_moho_top,ibelm_moho_bot,ibelm_400_top,ibelm_400_bot, &
@@ -747,6 +770,126 @@
!nspec_tiso = count(idoubling(1:nspec) == IFLAG_220_80) + count(idoubling(1:nspec) == IFLAG_80_MOHO)
nspec_tiso = count(ispec_is_tiso(:))
+!****************************************************************************************************
+! Mila
+
+ if(SORT_MESH_INNER_OUTER) then
+
+!!!! David Michea: detection of the edges, coloring and permutation separately
+ allocate(perm(nspec))
+
+! implement mesh coloring for GPUs if needed, to create subsets of disconnected elements
+! to remove dependencies and the need for atomic operations in the sum of elemental contributions in the solver
+ if(USE_MESH_COLORING_GPU) then
+
+ allocate(first_elem_number_in_this_color(MAX_NUMBER_OF_COLORS + 1))
+
+ call get_perm_color_faster(is_on_a_slice_edge,ibool,perm,nspec,nglob, &
+ nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,first_elem_number_in_this_color,myrank)
+
+! for the last color, the next color is fictitious and its first (fictitious) element number is nspec + 1
+ first_elem_number_in_this_color(nb_colors_outer_elements + nb_colors_inner_elements + 1) = nspec + 1
+
+ allocate(num_of_elems_in_this_color(nb_colors_outer_elements + nb_colors_inner_elements))
+
+! save mesh coloring
+ open(unit=99,file=prname(1:len_trim(prname))//'num_of_elems_in_this_color.dat',status='unknown')
+
+! number of colors for outer elements
+ write(99,*) nb_colors_outer_elements
+
+! number of colors for inner elements
+ write(99,*) nb_colors_inner_elements
+
+! number of elements in each color
+ do icolor = 1, nb_colors_outer_elements + nb_colors_inner_elements
+ num_of_elems_in_this_color(icolor) = first_elem_number_in_this_color(icolor+1) - first_elem_number_in_this_color(icolor)
+ write(99,*) num_of_elems_in_this_color(icolor)
+ enddo
+ close(99)
+
+! check that the sum of all the numbers of elements found in each color is equal
+! to the total number of elements in the mesh
+ if(sum(num_of_elems_in_this_color) /= nspec) then
+ print *,'nspec = ',nspec
+ print *,'total number of elements in all the colors of the mesh = ',sum(num_of_elems_in_this_color)
+ stop 'incorrect total number of elements in all the colors of the mesh'
+ endif
+
+! check that the sum of all the numbers of elements found in each color for the outer elements is equal
+! to the total number of outer elements found in the mesh
+ if(sum(num_of_elems_in_this_color(1:nb_colors_outer_elements)) /= nspec_outer) then
+ print *,'nspec_outer = ',nspec_outer
+ print *,'total number of elements in all the colors of the mesh for outer elements = ',sum(num_of_elems_in_this_color)
+ stop 'incorrect total number of elements in all the colors of the mesh for outer elements'
+ endif
+
+ call MPI_ALLREDUCE(nspec_outer,nspec_outer_min_global,1,MPI_INTEGER,MPI_MIN,MPI_COMM_WORLD,ier)
+ call MPI_ALLREDUCE(nspec_outer,nspec_outer_max_global,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,ier)
+
+ deallocate(first_elem_number_in_this_color)
+ deallocate(num_of_elems_in_this_color)
+
+ else
+
+!! DK DK for regular C + MPI version for CPUs: do not use colors but nonetheless put all the outer elements
+!! DK DK first in order to be able to overlap non-blocking MPI communications with calculations
+
+!! DK DK nov 2010, for Rosa Badia / StarSs:
+!! no need for mesh coloring, but need to implement inner/outer subsets for non blocking MPI for StarSs
+ ispec_counter = 0
+ perm(:) = 0
+
+! first generate all the outer elements
+ do ispec = 1,nspec
+ if(is_on_a_slice_edge(ispec)) then
+ ispec_counter = ispec_counter + 1
+ perm(ispec) = ispec_counter
+ endif
+ enddo
+
+! make sure we have detected some outer elements
+ if(ispec_counter <= 0) stop 'fatal error: no outer elements detected!'
+
+! store total number of outer elements
+ nspec_outer = ispec_counter
+
+! then generate all the inner elements
+ do ispec = 1,nspec
+ if(.not. is_on_a_slice_edge(ispec)) then
+ ispec_counter = ispec_counter + 1
+ perm(ispec) = ispec_counter
+ endif
+ enddo
+
+! test that all the elements have been used once and only once
+ if(ispec_counter /= nspec) stop 'fatal error: ispec_counter not equal to nspec'
+
+! do basic checks
+ if(minval(perm) /= 1) stop 'minval(perm) should be 1'
+ if(maxval(perm) /= nspec) stop 'maxval(perm) should be nspec'
+
+ call MPI_ALLREDUCE(nspec_outer,nspec_outer_min_global,1,MPI_INTEGER,MPI_MIN,MPI_COMM_WORLD,ier)
+ call MPI_ALLREDUCE(nspec_outer,nspec_outer_max_global,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,ier)
+
+ endif
+
+ else
+!
+ print *,'SORT_MESH_INNER_OUTER must always been set to .true. even for the regular C version for CPUs'
+ print *,'in order to be able to use non blocking MPI to overlap communications'
+! print *,'generating identity permutation'
+! do ispec = 1,nspec
+! perm(ispec) = ispec
+! enddo
+ stop 'please set SORT_MESH_INNER_OUTER to .true. and recompile the whole code'
+
+ endif
+
+!!!! David Michea: end of mesh coloring
+
+!****************************************************************************************************
+
! precomputes jacobian for 2d absorbing boundary surfaces
call get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90 2011-06-15 22:59:57 UTC (rev 18639)
@@ -0,0 +1,813 @@
+
+! define sets of colors that contain disconnected elements for the CUDA solver.
+! also split the elements into two subsets: inner and outer elements, in order
+! to be able to compute the outer elements first in the solver and then
+! start non-blocking MPI calls and overlap them with the calculation of the inner elements
+! (which works fine because there are always far more inner elements than outer elements)
+
+!*********************************************************************************************************
+! Mila
+subroutine get_perm_color_faster(is_on_a_slice_edge,ibool,perm,nspec,nglob, &
+ nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,first_elem_number_in_this_color,myrank)
+
+ implicit none
+
+ include "constants.h"
+
+ logical, dimension(nspec) :: is_on_a_slice_edge
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ integer, dimension(nspec) :: perm
+ integer, dimension(nspec) :: color
+ integer, dimension(MAX_NUMBER_OF_COLORS) :: first_elem_number_in_this_color
+ integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,myrank
+
+! local variables
+ integer nspec, nglob
+
+ call get_color_faster(ibool, is_on_a_slice_edge, myrank, nspec, nglob, &
+ color, nb_colors_outer_elements, nb_colors_inner_elements, nspec_outer)
+
+ if(myrank == 0) then
+ write(IMAIN,*) 'number of colors of the graph for inner elements = ',nb_colors_inner_elements
+ write(IMAIN,*) 'number of colors of the graph for outer elements = ',nb_colors_outer_elements
+ write(IMAIN,*) 'total number of colors of the graph (sum of both) = ', nb_colors_inner_elements + nb_colors_outer_elements
+ write(IMAIN,*) 'number of elements of the graph for outer elements = ',nspec_outer
+ endif
+
+ if(myrank == 0) write(IMAIN,*) 'generating the final colors'
+ first_elem_number_in_this_color(:) = -1
+ call get_final_perm(color,perm,first_elem_number_in_this_color,nspec,nb_colors_inner_elements+nb_colors_outer_elements)
+
+ if(myrank == 0) write(IMAIN,*) 'done with mesh coloring and inner/outer element splitting'
+
+ end subroutine get_perm_color_faster
+
+!-----------------------------------------------------------------------
+
+subroutine get_color_faster(ibool, is_on_a_slice_edge, myrank, nspec, nglob, &
+ color, nb_colors_outer_elements, nb_colors_inner_elements, nspec_outer)
+
+ implicit none
+
+ include "constants.h"
+
+ logical, dimension(nspec) :: is_on_a_slice_edge
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ integer, dimension(nspec) :: color
+ integer :: nb_colors_outer_elements,nb_colors_inner_elements,myrank,nspec_outer
+
+! local variables
+ integer nspec,nglob
+
+ integer :: ispec
+
+!! DK DK for mesh coloring GPU Joseph Fourier
+ logical, dimension(:), allocatable :: mask_ibool
+
+ integer :: icolor, nb_already_colored
+ integer :: iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
+
+ logical :: conflict_found_need_new_color
+
+ nspec_outer = 0
+ do ispec=1,nspec
+ if(is_on_a_slice_edge(ispec)) nspec_outer=nspec_outer+1
+ enddo
+
+!! DK DK start mesh coloring (new Apr 2010 version by DK for GPU Joseph Fourier)
+ allocate(mask_ibool(nglob))
+
+!! DK DK ----------------------------------
+!! DK DK color the mesh in the crust_mantle
+!! DK DK ----------------------------------
+
+ if(myrank == 0) then
+ print *
+ print *,'----------------------------------'
+ print *,'coloring the mesh'
+ print *,'----------------------------------'
+ print *
+ endif
+
+
+! first set color of all elements to 1
+ color(:) = 0
+ icolor = 0
+ nb_already_colored = 0
+
+ do while(nb_already_colored<nspec_outer)
+ icolor = icolor + 1
+ 333 continue
+ !if(myrank == 0) print *,'analyzing color ',icolor,' for all the elements of the mesh'
+ mask_ibool(:) = .false.
+ conflict_found_need_new_color = .false.
+
+ do ispec = 1,nspec
+ if(is_on_a_slice_edge(ispec)) then
+ if(color(ispec) == 0) then
+ ! the eight corners of the current element
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(mask_ibool(iglob1) .or. mask_ibool(iglob2) .or. mask_ibool(iglob3) .or. mask_ibool(iglob4) .or. &
+ mask_ibool(iglob5) .or. mask_ibool(iglob6) .or. mask_ibool(iglob7) .or. mask_ibool(iglob8)) then
+ ! if element of this color has a common point with another element of that same color
+ ! then we need to create a new color, i.e., increment the color of the current element
+ conflict_found_need_new_color = .true.
+ else
+ color(ispec) = icolor
+ nb_already_colored = nb_already_colored + 1
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob5) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ endif
+ endif
+ endif
+ enddo
+
+ if(conflict_found_need_new_color) then
+ icolor = icolor + 1
+ goto 333
+ endif
+ enddo
+
+ nb_colors_outer_elements = icolor
+
+ do while(nb_already_colored<nspec)
+ icolor = icolor + 1
+ 334 continue
+ !if(myrank == 0) print *,'analyzing color ',icolor,' for all the elements of the mesh'
+ mask_ibool(:) = .false.
+ conflict_found_need_new_color = .false.
+ do ispec = 1,nspec
+ if (.not. is_on_a_slice_edge(ispec)) then
+ if(color(ispec) == 0) then
+ ! the eight corners of the current element
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(mask_ibool(iglob1) .or. mask_ibool(iglob2) .or. mask_ibool(iglob3) .or. mask_ibool(iglob4) .or. &
+ mask_ibool(iglob5) .or. mask_ibool(iglob6) .or. mask_ibool(iglob7) .or. mask_ibool(iglob8)) then
+ ! if element of this color has a common point with another element of that same color
+ ! then we need to create a new color, i.e., increment the color of the current element
+ conflict_found_need_new_color = .true.
+ else
+ color(ispec) = icolor
+ nb_already_colored = nb_already_colored + 1
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob5) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ endif
+ endif
+ endif
+ enddo
+
+ if(conflict_found_need_new_color) then
+ icolor = icolor + 1
+ goto 334
+ endif
+ enddo
+
+ nb_colors_inner_elements = icolor - nb_colors_outer_elements
+
+ if(myrank == 0) print *,'created a total of ',maxval(color),' colors for all the elements of the mesh'
+
+ !!!!!!!! DK DK now check that all the color sets are independent
+ do icolor = 1,maxval(color)
+ mask_ibool(:) = .false.
+ do ispec = 1,nspec
+ if(color(ispec) == icolor) then
+ ! the eight corners of the current element
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(mask_ibool(iglob1) .or. mask_ibool(iglob2) .or. mask_ibool(iglob3) .or. mask_ibool(iglob4) .or. &
+ mask_ibool(iglob5) .or. mask_ibool(iglob6) .or. mask_ibool(iglob7) .or. mask_ibool(iglob8)) then
+ ! if element of this color has a common point with another element of that same color
+ ! then there is a problem, the color set is not correct
+ stop 'error detected: found a common point inside a color set'
+ else
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob5) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ endif
+ endif
+ enddo
+
+ if(myrank == 0) print *,'color ',icolor,' has disjoint elements only and is therefore OK'
+ if(myrank == 0) print *,'it contains ',count(color == icolor),' elements'
+ enddo
+
+ if(myrank == 0) print *,'the ',maxval(color),' color sets are OK'
+
+ deallocate(mask_ibool)
+
+end subroutine get_color_faster
+
+!*********************************************************************************************************
+
+ subroutine get_perm_color(is_on_a_slice_edge,ibool,perm,nspec,nglob, &
+ nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,first_elem_number_in_this_color,myrank)
+
+ implicit none
+
+ include "constants.h"
+
+ logical, dimension(nspec) :: is_on_a_slice_edge
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ integer, dimension(nspec) :: perm
+ integer, dimension(nspec) :: color
+ integer, dimension(MAX_NUMBER_OF_COLORS) :: first_elem_number_in_this_color
+ integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,myrank
+
+! local variables
+ integer nspec,nglob_GLL_full
+
+! a neighbor of a hexahedral node is a hexahedron that shares a face with it -> max degree of a node = 6
+ integer, parameter :: MAX_NUMBER_OF_NEIGHBORS = 100
+
+! global corner numbers that need to be created
+ integer, dimension(nglob) :: global_corner_number
+
+ integer mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
+ integer, dimension(:), allocatable :: ne,np,adj
+ integer xadj(nspec+1)
+
+ !logical maskel(nspec)
+
+ integer i,istart,istop,number_of_neighbors
+
+ integer nglob_eight_corners_only,nglob
+
+! only count the total size of the array that will be created, or actually create it
+ logical count_only
+ integer total_size_ne,total_size_adj
+
+!
+!-----------------------------------------------------------------------
+!
+
+! total number of points in the mesh
+ nglob_GLL_full = nglob
+
+!---- call Charbel Farhat's routines
+ if(myrank == 0) &
+ write(IMAIN,*) 'calling form_elt_connectivity_foelco to perform mesh coloring and inner/outer element splitting'
+ call form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number,nglob_GLL_full,ibool,nglob_eight_corners_only)
+ do i=1,nspec
+ istart = mp(i)
+ istop = mp(i+1) - 1
+ enddo
+
+! count only, to determine the size needed for the array
+ allocate(np(nglob_eight_corners_only+1))
+ count_only = .true.
+ total_size_ne = 1
+ if(myrank == 0) write(IMAIN,*) 'calling form_node_connectivity_fonoco to determine the size of the table'
+ allocate(ne(total_size_ne))
+ call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
+ deallocate(ne)
+
+!print *, 'nglob_eight_corners_only'
+!print *, nglob_eight_corners_only
+
+! allocate the array with the right size
+ allocate(ne(total_size_ne))
+
+! now actually generate the array
+ count_only = .false.
+ if(myrank == 0) write(IMAIN,*) 'calling form_node_connectivity_fonoco to actually create the table'
+ call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
+ do i=1,nglob_eight_corners_only
+ istart = np(i)
+ istop = np(i+1) - 1
+ enddo
+
+!print *, 'total_size_ne'
+!print *, total_size_ne
+
+! count only, to determine the size needed for the array
+ count_only = .true.
+ total_size_adj = 1
+ if(myrank == 0) write(IMAIN,*) 'calling create_adjacency_table_adjncy to determine the size of the table'
+ allocate(adj(total_size_adj))
+ !call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
+ !count_only,total_size_ne,total_size_adj,.false.)
+ call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,nspec,nglob_eight_corners_only,&
+ count_only,total_size_ne,total_size_adj,.false.)
+ deallocate(adj)
+
+! allocate the array with the right size
+ allocate(adj(total_size_adj))
+
+! now actually generate the array
+ count_only = .false.
+ if(myrank == 0) write(IMAIN,*) 'calling create_adjacency_table_adjncy again to actually create the table'
+ !call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
+ !count_only,total_size_ne,total_size_adj,.false.)
+ call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,nspec,nglob_eight_corners_only,&
+ count_only,total_size_ne,total_size_adj,.false.)
+
+ do i=1,nspec
+ istart = xadj(i)
+ istop = xadj(i+1) - 1
+ number_of_neighbors = istop-istart+1
+ if(number_of_neighbors < 1 .or. number_of_neighbors > MAX_NUMBER_OF_NEIGHBORS) stop 'incorrect number of neighbors'
+ enddo
+
+ deallocate(ne,np)
+
+ call get_color(adj,xadj,color,nspec,total_size_adj,is_on_a_slice_edge, &
+ nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer)
+
+ if(myrank == 0) then
+ write(IMAIN,*) 'number of colors of the graph for inner elements = ',nb_colors_inner_elements
+ write(IMAIN,*) 'number of colors of the graph for outer elements = ',nb_colors_outer_elements
+ write(IMAIN,*) 'total number of colors of the graph (sum of both) = ', nb_colors_inner_elements + nb_colors_outer_elements
+ write(IMAIN,*) 'number of elements of the graph for outer elements = ',nspec_outer
+ endif
+
+ deallocate(adj)
+
+ if(myrank == 0) write(IMAIN,*) 'generating the final colors'
+ first_elem_number_in_this_color(:) = -1
+ call get_final_perm(color,perm,first_elem_number_in_this_color,nspec,nb_colors_inner_elements+nb_colors_outer_elements)
+
+ if(myrank == 0) write(IMAIN,*) 'done with mesh coloring and inner/outer element splitting'
+
+ end subroutine get_perm_color
+
+!------------------------------------------------------------------
+
+subroutine get_final_perm(color,perm,first_elem_number_in_this_color,nspec,nb_color)
+
+ integer, intent(in) :: nspec,nb_color
+ integer, intent(in) :: color(nspec)
+ integer, intent(inout) :: perm(nspec)
+ integer, intent(inout) :: first_elem_number_in_this_color(nb_color)
+ integer :: ielem,icolor,counter
+
+ counter = 1
+ do icolor = 1, nb_color
+ first_elem_number_in_this_color(icolor) = counter
+ do ielem = 1, nspec
+ if(color(ielem) == icolor) then
+ perm(ielem) = counter
+ counter = counter + 1
+ endif
+ enddo
+ enddo
+
+end subroutine get_final_perm
+
+!------------------------------------------------------------------
+
+subroutine get_color(adj,xadj,color,nspec,total_size_adj,is_on_a_slice_edge, &
+ nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer)
+
+ integer, intent(in) :: nspec,total_size_adj
+ integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
+ integer :: color(nspec)
+ integer :: this_color,nb_already_colored,ispec,ixadj,ok
+ logical, dimension(nspec) :: is_on_a_slice_edge
+ integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer
+ logical :: is_outer_element(nspec)
+
+ nspec_outer = 0
+
+ is_outer_element(:) = .false.
+
+ do ispec=1,nspec
+ if (is_on_a_slice_edge(ispec)) then
+ is_outer_element(ispec) = .true.
+ nspec_outer=nspec_outer+1
+ endif
+ enddo
+
+! outer elements
+ color(:) = 0
+ this_color = 0
+ nb_already_colored = 0
+ do while(nb_already_colored<nspec_outer)
+ this_color = this_color + 1
+ do ispec = 1, nspec
+ if (is_outer_element(ispec)) then
+ if (color(ispec) == 0) then
+ ok = 1
+ do ixadj = xadj(ispec), (xadj(ispec+1)-1)
+ if (is_outer_element(adj(ixadj)) .and. color(adj(ixadj)) == this_color) ok = 0
+ enddo
+ if (ok /= 0) then
+ color(ispec) = this_color
+ nb_already_colored = nb_already_colored + 1
+ endif
+ endif
+ endif
+ enddo
+ enddo
+ nb_colors_outer_elements = this_color
+
+! inner elements
+ do while(nb_already_colored<nspec)
+ this_color = this_color + 1
+ do ispec = 1, nspec
+ if (.not. is_outer_element(ispec)) then
+ if (color(ispec) == 0) then
+ ok = 1
+ do ixadj = xadj(ispec), (xadj(ispec+1)-1)
+ if (.not. is_outer_element(adj(ixadj)) .and. color(adj(ixadj)) == this_color) ok = 0
+ enddo
+ if (ok /= 0) then
+ color(ispec) = this_color
+ nb_already_colored = nb_already_colored + 1
+ endif
+ endif
+ endif
+ enddo
+ enddo
+
+ nb_colors_inner_elements = this_color - nb_colors_outer_elements
+
+end subroutine get_color
+
+!------------------------------------------------------------------
+
+!=======================================================================
+!
+! Charbel Farhat's FEM topology routines
+!
+! Dimitri Komatitsch, February 1996 - Code based on Farhat's original version from 1987
+!
+! modified and adapted by Dimitri Komatitsch, May 2006
+!
+!=======================================================================
+
+ subroutine form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number,&
+nglob_GLL_full,ibool,nglob_eight_corners_only)
+
+!-----------------------------------------------------------------------
+!
+! Forms the MN and MP arrays
+!
+! Input :
+! -------
+! ibool Array needed to build the element connectivity table
+! nspec Number of elements in the domain
+! NGNOD_HEXAHEDRA number of nodes per hexahedron (brick with 8 corners)
+!
+! Output :
+! --------
+! MN, MP This is the element connectivity array pair.
+! Array MN contains the list of the element
+! connectivity, that is, the nodes contained in each
+! element. They are stored in a stacked fashion.
+!
+! Pointer array MP stores the location of each
+! element list. Its length is equal to the number
+! of elements plus one.
+!
+!-----------------------------------------------------------------------
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,nglob_GLL_full
+
+! arrays with mesh parameters per slice
+ integer, intent(in), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! global corner numbers that need to be created
+ integer, intent(out), dimension(nglob_GLL_full) :: global_corner_number
+ integer, intent(out) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
+ integer, intent(out) :: nglob_eight_corners_only
+
+ integer ninter,nsum,ispec,node,k,inumcorner,ix,iy,iz
+
+ ninter = 1
+ nsum = 1
+ mp(1) = 1
+
+!---- define topology of the elements in the mesh
+!---- we need to define adjacent numbers from the sub-mesh consisting of the corners only
+ nglob_eight_corners_only = 0
+ global_corner_number(:) = -1
+
+ do ispec=1,nspec
+
+ inumcorner = 0
+ do iz = 1,NGLLZ,NGLLZ-1
+ do iy = 1,NGLLY,NGLLY-1
+ do ix = 1,NGLLX,NGLLX-1
+
+ inumcorner = inumcorner + 1
+ if(inumcorner > NGNOD_HEXAHEDRA) stop 'corner number too large'
+
+! check if this point was already assigned a number previously, otherwise create one and store it
+ if(global_corner_number(ibool(ix,iy,iz,ispec)) == -1) then
+ nglob_eight_corners_only = nglob_eight_corners_only + 1
+ global_corner_number(ibool(ix,iy,iz,ispec)) = nglob_eight_corners_only
+ endif
+
+ node = global_corner_number(ibool(ix,iy,iz,ispec))
+ do k=nsum,ninter-1
+ if(node == mn(k)) goto 200
+ enddo
+
+ mn(ninter) = node
+ ninter = ninter + 1
+ 200 continue
+
+ enddo
+ enddo
+ enddo
+
+ nsum = ninter
+ mp(ispec + 1) = nsum
+
+ enddo
+
+ end subroutine form_elt_connectivity_foelco
+
+!
+!----------------------------------------------------
+!
+
+ subroutine form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,&
+nspec,count_only,total_size_ne)
+
+!-----------------------------------------------------------------------
+!
+! Forms the NE and NP arrays
+!
+! Input :
+! -------
+! MN, MP, nspec
+! nglob_eight_corners_only Number of nodes in the domain
+!
+! Output :
+! --------
+! NE, NP This is the node-connected element array pair.
+! Integer array NE contains a list of the
+! elements connected to each node, stored in stacked fashion.
+!
+! Array NP is the pointer array for the
+! location of a node's element list in the NE array.
+! Its length is equal to the number of points plus one.
+!
+!-----------------------------------------------------------------------
+
+ implicit none
+
+ include "constants.h"
+
+! only count the total size of the array that will be created, or actually create it
+ logical count_only
+ integer total_size_ne
+
+ integer nglob_eight_corners_only,nspec
+
+ integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
+
+ integer, intent(out) :: ne(total_size_ne),np(nglob_eight_corners_only+1)
+
+ integer nsum,inode,ispec,j
+
+ nsum = 1
+ np(1) = 1
+
+ do inode=1,nglob_eight_corners_only
+ do 200 ispec=1,nspec
+
+ do j=mp(ispec),mp(ispec + 1) - 1
+ if (mn(j) == inode) then
+ if(count_only) then
+ total_size_ne = nsum
+ else
+ ne(nsum) = ispec
+ endif
+ nsum = nsum + 1
+ goto 200
+ endif
+ enddo
+ 200 continue
+
+ np(inode + 1) = nsum
+
+ enddo
+
+ end subroutine form_node_connectivity_fonoco
+
+!
+!----------------------------------------------------
+!
+
+ !subroutine create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
+ ! count_only,total_size_ne,total_size_adj,face)
+ subroutine create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,nspec,nglob_eight_corners_only,&
+ count_only,total_size_ne,total_size_adj,face)
+
+!-----------------------------------------------------------------------
+!
+! Establishes the element adjacency information of the mesh
+! Two elements are considered adjacent if they share a face.
+!
+! Input :
+! -------
+! MN, MP, NE, NP, nspec
+! MASKEL logical mask (length = nspec)
+!
+! Output :
+! --------
+! ADJ, XADJ This is the element adjacency array pair. Array
+! ADJ contains the list of the elements adjacent to
+! element i. They are stored in a stacked fashion.
+! Pointer array XADJ stores the location of each element list.
+!
+!-----------------------------------------------------------------------
+
+ implicit none
+
+ include "constants.h"
+
+! only count the total size of the array that will be created, or actually create it
+ logical count_only,face
+ integer total_size_ne,total_size_adj
+
+ integer nglob_eight_corners_only
+
+ integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1),ne(total_size_ne),np(nglob_eight_corners_only+1)
+
+ integer, intent(out) :: adj(total_size_adj),xadj(nspec+1)
+
+ logical maskel(nspec)
+ integer countel(nspec)
+
+ integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel
+
+ xadj(1) = 1
+ iad = 1
+
+ do ispec=1,nspec
+
+! reset mask
+ maskel(:) = .false.
+
+! mask current element
+ maskel(ispec) = .true.
+ if (face) countel(:) = 0
+
+ istart = mp(ispec)
+ istop = mp(ispec+1) - 1
+ do ino=istart,istop
+ node = mn(ino)
+ jstart = np(node)
+ jstop = np(node + 1) - 1
+ do 120 jel=jstart,jstop
+ nelem = ne(jel)
+ if(maskel(nelem)) goto 120
+ if (face) then
+ ! if 2 elements share at least 3 corners, therefore they share a face
+ countel(nelem) = countel(nelem) + 1
+ if (countel(nelem)>=3) then
+ if(count_only) then
+ total_size_adj = iad
+ else
+ adj(iad) = nelem
+ endif
+ maskel(nelem) = .true.
+ iad = iad + 1
+ endif
+ else
+ if(count_only) then
+ total_size_adj = iad
+ else
+ adj(iad) = nelem
+ endif
+ maskel(nelem) = .true.
+ iad = iad + 1
+ endif
+ 120 continue
+ enddo
+
+ xadj(ispec+1) = iad
+
+ enddo
+
+ end subroutine create_adjacency_table_adjncy
+
+
+ subroutine permute_elements_real(array_to_permute,temp_array,perm,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: nspec
+ integer, intent(in), dimension(nspec) :: perm
+
+ real(kind=CUSTOM_REAL), intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
+
+ integer old_ispec,new_ispec
+
+! copy the original array
+ temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
+
+ do old_ispec = 1,nspec
+ new_ispec = perm(old_ispec)
+ array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
+ enddo
+
+ end subroutine permute_elements_real
+
+!
+!-----------------------------------------------------------------------
+!
+
+! implement permutation of elements for arrays of integer type
+ subroutine permute_elements_integer(array_to_permute,temp_array,perm,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: nspec
+ integer, intent(in), dimension(nspec) :: perm
+
+ integer, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
+
+ integer old_ispec,new_ispec
+
+! copy the original array
+ temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
+
+ do old_ispec = 1,nspec
+ new_ispec = perm(old_ispec)
+ array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
+ enddo
+
+ end subroutine permute_elements_integer
+
+!
+!-----------------------------------------------------------------------
+!
+
+! implement permutation of elements for arrays of double precision type
+ subroutine permute_elements_dble(array_to_permute,temp_array,perm,nspec)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: nspec
+ integer, intent(in), dimension(nspec) :: perm
+
+ double precision, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
+
+ integer old_ispec,new_ispec
+
+! copy the original array
+ temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
+
+ do old_ispec = 1,nspec
+ new_ispec = perm(old_ispec)
+ array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
+ enddo
+
+ end subroutine permute_elements_dble
More information about the CIG-COMMITS
mailing list