[cig-commits] r18320 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Thu May 5 00:38:54 PDT 2011
Author: percygalvez
Date: 2011-05-05 00:38:54 -0700 (Thu, 05 May 2011)
New Revision: 18320
Removed:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/get_coupling_domain1_domain2.f90
Log:
deleting obsolete files
Deleted: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/get_coupling_domain1_domain2.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/get_coupling_domain1_domain2.f90 2011-05-05 05:20:33 UTC (rev 18319)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/get_coupling_domain1_domain2.f90 2011-05-05 07:38:54 UTC (rev 18320)
@@ -1,408 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-module coupling
-
- implicit none
- type (coupling_type)
- private
- double precision, dimension(:,:) :: jacobian2Dw
- double precision, dimension(:,:,:) :: normal,ijk
-
- integer :: tag1,tag2,num_faces,NbFaults
- integer, dimension(:) :: ispec
- end type coupling_type
- type (coupling_type),pointer :: fault_db(:)
- public :: get_coupling_surfaces_domain2_domain1
-
-contains
-
-!=====================================================================
-subroutine read_parameters_fault
-
- implicit none
-
-! open Par_file_fault
- open(unit=100,file='~/SPECFEM3D_FAULT/DATA/Par_file_faults.in')
-! if file does not exist: NbFaults=0
-
- read(*,*) fault_db%NbFaults
-! if already allocated (associated), deallocate
- allocate(fault_db%NbFaults)
- do i=1,fault_db%NbFaults
- read(*,*) fault_db(i)%tag1,fault_db(i)%tag2
- enddo
-
-! close Par_file_fault
- close(100)
-
-
-end subroutine read_parameters_fault
-
-!=====================================================================
-
-
-
-!=====================================================================
-! BEGIN INPUT BLOCK
-
-! inputs:
-! myrank: processor index
-! domain_1tag and domain_2tag : tags created by cubit specifying different domain.
-! nspec : number of spectral elements in each block.
-! nglob : number of gobal nodes in each block.
-! ibool : local to global numbering table, iglob = ibool(i,j,k,ispec)
-
-! do inum = 1,coupling%num_faces
-! coupling%normal(:,:,inum) = tmp_normal(:,:,inum)
-! coupling%jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
-! coupling%ijk(:,:,inum) = tmp_ijk(:,:,inum)
-! coupling%ispec(inum) = tmp_ispec(inum)
-!
-! INPUTS for MPI comunications
-! NPROC : number of processors.
-! nibool_interfaces_ext_mesh:
-! ibool_interfaces_ext_mesh:
-! num_interfaces_ext_mesh :
-! max_interface_size_ext_mesh:
-! my_neighbours_ext_mesh:
-
-!
-! OUTPUTS:
-! coupling: fault structure (database)
-
-
-! determines coupling surface for domain2-domain1 domains
-
-subroutine get_coupling_surfaces_domain2_domain1(myrank, &
- coupling,domain1_flag,domain2_flag, &
- nspec,nglob,ibool,NPROC, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
- num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
- my_neighbours_ext_mesh)
-
-
- use generate_databases_par, only:mat_ext_mesh,elmnts_ext_mesh ! mat_ext_mesh , elemnts_ext_mesh
- use create_regions_mesh_ext_par
-
- type(coupling_type), intent(inout) ::coupling
- integer, dimension(:), allocatable :: domain1_flag,domain2_flag
-
-
-! domain1_flag=fault_db%tag1
-! domain2_flag=fault_db%tag2
-
-
-! number of spectral elements in each block
- integer :: myrank,nspec,nglob,NPROC
-
-! arrays with the mesh
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! MPI communication
- integer :: num_interfaces_ext_mesh,max_interface_size_ext_mesh
- integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
- integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: &
- ibool_interfaces_ext_mesh
- integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
-
-! local parameters
- ! (assumes NGLLX=NGLLY=NGLLZ)
-! only defining local variables...
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
- real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
- real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
- real(kind=CUSTOM_REAL),dimension(:,:,:),allocatable :: tmp_normal
- real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: tmp_jacobian2Dw
- integer :: ijk_face(3,NGLLX,NGLLY)
- integer,dimension(:,:,:),allocatable :: tmp_ijk
- integer,dimension(:),allocatable :: tmp_ispec
- integer,dimension(NGNOD2D) :: iglob_corners_ref !,iglob_corners
- integer :: ispec,i,j,k,igll,ier,iglob
- integer :: inum,iface_ref,icorner,iglob_midpoint ! iface,ispec_neighbor
- integer :: count_domain1,count_domain2
- ! mpi interface communication
- integer, dimension(:), allocatable :: test_flag
- integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
- integer :: max_nibool_interfaces_ext_mesh
- logical, dimension(:), allocatable :: mask_ibool
-
- ! corners indices of reference cube faces
- integer,dimension(3,4),parameter :: iface1_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
- integer,dimension(3,4),parameter :: iface2_corner_ijk = &
- reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
- integer,dimension(3,4),parameter :: iface3_corner_ijk = &
- reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
- integer,dimension(3,4),parameter :: iface4_corner_ijk = &
- reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
- integer,dimension(3,4),parameter :: iface5_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
- integer,dimension(3,4),parameter :: iface6_corner_ijk = &
- reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
- integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
- reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
- iface3_corner_ijk,iface4_corner_ijk, &
- iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
- integer,dimension(3,6),parameter :: iface_all_midpointijk = &
- reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
-
-
- ! test vtk output
- !integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
- !character(len=256):: prname_file
-
-
-! allocates temporary arrays
- allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6))
- allocate(tmp_jacobian2Dw(NGLLSQUARE,nspec*6))
- allocate(tmp_ijk(3,NGLLSQUARE,nspec*6))
- allocate(tmp_ispec(nspec*6))
- tmp_ispec(:) = 0
- tmp_ijk(:,:,:) = 0
- tmp_normal(:,:,:) = 0.0
- tmp_jacobian2Dw(:,:) = 0.0
-
- ! sets flags for domain2 / domain1 on global points
- allocate(domain1_flag(nglob),stat=ier)
- allocate(domain2_flag(nglob),stat=ier)
-
- ! what is this test_flag is for .
- allocate(test_flag(nglob),stat=ier)
- allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocate flag array'
- domain1_flag(:) = 0
- domain2_flag(:) = 0
- test_flag(:) = 0
- count_domain1 = 0
- count_domain2 = 0
-
-!!!! running onver all the elements over each block.
-
- do ispec = 1, nspec
- ! counts elements
-!!!change this variable.. my_tag(ispec)==domain2+_tag)
-!!!change this variable.. my_tag(ispec)==domain2-_tag)
-!!!! allocate my_tag
-! my_tag = mat_ext_mesh(1,ispec)
-
- my_tag(ispec)=mat_ext_mesh(1,ispec)
- if( my_tag(ispec)==domain1_tag ) count_domain1 = count_domain1 + 1
- if( my_tag(ispec)==domain2_tag ) count_domain2 = count_domain2 + 1
-
-!!!!!! inserting domains into processor , one by one. ...
- ! sets flags on global points
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- ! global index
- iglob = ibool(i,j,k,ispec)
- ! sets domain1 flag1
- if( ispec_is_domain1(ispec) ) domain1_flag(iglob) = myrank+1
- ! sets domain2 flag2
- if( ispec_is_domain2(ispec) ) domain2_flag(iglob) = myrank+1
- ! sets test flag
- test_flag(iglob) = myrank+1
- enddo
- enddo
- enddo
- enddo
-
-!!! counting number of domains and assigning them into each processors.
- call sum_all_i(count_domain2,inum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' total domain2 elements:',inum
- endif
- call sum_all_i(count_domain1,inum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' total domain1 elements :',inum
- endif
-
-
-
- ! collects contributions from different MPI partitions
- ! sets up MPI communications
- max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
- allocate(ibool_interfaces_ext_mesh_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
- if( ier /= 0 ) stop 'error allocating array'
- do i = 1, num_interfaces_ext_mesh
- ibool_interfaces_ext_mesh_dummy(:,i) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,i)
- enddo
- ! sums domain1 flags
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,domain1_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
- my_neighbours_ext_mesh)
- ! sums domain2 flags
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,domain2_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
- my_neighbours_ext_mesh)
-
- ! sums test flags
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,test_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
- my_neighbours_ext_mesh)
-
- ! loops over all element faces and
- ! counts number of coupling faces between domain2 and domain1 elements
- mask_ibool(:) = .false.
- inum = 0
- do ispec=1,nspec
-
- ! loops over each face
- do iface_ref= 1, 6
-
- ! takes indices of corners of reference face
- do icorner = 1,NGNOD2D
- i = iface_all_corner_ijk(1,icorner,iface_ref)
- j = iface_all_corner_ijk(2,icorner,iface_ref)
- k = iface_all_corner_ijk(3,icorner,iface_ref)
- ! global reference indices
- iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
-
- ! reference corner coordinates
- xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
- ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
- zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
- enddo
-
- ! checks if face has domain2 side
- if( domain2_flag( iglob_corners_ref(1) ) >= 1 .and. &
- domain2_flag( iglob_corners_ref(2) ) >= 1 .and. &
- domain2_flag( iglob_corners_ref(3) ) >= 1 .and. &
- domain2_flag( iglob_corners_ref(4) ) >= 1) then
- ! checks if face is has an domain1 side
- if( domain1_flag( iglob_corners_ref(1) ) >= 1 .and. &
- domain1_flag( iglob_corners_ref(2) ) >= 1 .and. &
- domain1_flag( iglob_corners_ref(3) ) >= 1 .and. &
- domain1_flag( iglob_corners_ref(4) ) >= 1) then
-
- ! reference midpoint on face (used to avoid redundant face counting)
- i = iface_all_midpointijk(1,iface_ref)
- j = iface_all_midpointijk(2,iface_ref)
- k = iface_all_midpointijk(3,iface_ref)
- iglob_midpoint = ibool(i,j,k,ispec)
-
- ! checks if points on this face are masked already
- if( .not. mask_ibool(iglob_midpoint) ) then
-
- ! gets face GLL points i,j,k indices from element face
- call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
-
- ! takes each element face only once, if it lies on an MPI interface
- ! note: this is not exactly load balanced
- ! lowest rank process collects as many faces as possible, second lowest as so forth
- if( (test_flag(iglob_midpoint) == myrank+1) .or. &
- (test_flag(iglob_midpoint) > 2*(myrank+1)) ) then
-
- ! gets face GLL 2Djacobian, weighted from element face
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
-
- ! normal convention: points away from domain2, reference element
- ! switch normal direction if necessary
- do j=1,NGLLY
- do i=1,NGLLX
- ! directs normals such that they point outwards of element
- call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- ! makes sure that it always points away from domain2 element,
- ! otherwise switch direction
- if( ispec_is_domain1(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
- enddo
- enddo
-
- ! stores informations about this face
- inum = inum + 1
- tmp_ispec(inum) = ispec
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- ! adds all gll points on this face
- igll = igll + 1
-
- ! do we need to store local i,j,k,ispec info? or only global indices iglob?
- tmp_ijk(:,igll,inum) = ijk_face(:,i,j)
-
- ! stores weighted jacobian and normals
- tmp_jacobian2Dw(igll,inum) = jacobian2Dw_face(i,j)
- tmp_normal(:,igll,inum) = normal_face(:,i,j)
-
- ! masks global points ( to avoid redundant counting of faces)
- iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
- mask_ibool(iglob) = .true.
- enddo
- enddo
- else
- ! assumes to be already collected by lower rank process, masks face points
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
- mask_ibool(iglob) = .true.
- enddo
- enddo
- endif ! test_flag
- endif ! mask_ibool
- endif ! domain1_flag
- endif ! domain2_flag
- enddo ! iface_ref
- enddo ! ispec
-
-! stores completed coupling domain2-domain1 face informations
-!
-! note: no need to store material parameters on these coupling points
-! for domain2-domain1 interface
-! defining new parameter , renaming coupling_fa_el_normal ..
-
-
- coupling%num_faces = inum
- allocate(coupling%normal(NDIM,NGLLSQUARE,coupling%num_faces))
- allocate(coupling%jacobian2Dw(NGLLSQUARE,coupling%num_faces))
- allocate(coupling%ijk(3,NGLLSQUARE,coupling%num_faces))
- allocate(coupling%ispec(coupling%num_faces))
- do inum = 1,coupling%num_faces
- coupling%normal(:,:,inum) = tmp_normal(:,:,inum)
- coupling%jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
- coupling%ijk(:,:,inum) = tmp_ijk(:,:,inum)
- coupling%ispec(inum) = tmp_ispec(inum)
- enddo
-
-! user output
- call sum_all_i(coupling%num_faces,inum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' domain2-domain1 coupling:'
- write(IMAIN,*) ' total number of faces = ',inum
- endif
-
- end subroutine get_coupling_surfaces_domain2_domain1
-
More information about the CIG-COMMITS
mailing list