[cig-commits] r18366 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: decompose_mesh_SCOTCH src
ampuero at geodynamics.org
ampuero at geodynamics.org
Fri May 13 20:08:08 PDT 2011
Author: ampuero
Date: 2011-05-13 20:08:08 -0700 (Fri, 13 May 2011)
New Revision: 18366
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
Log:
cleaned up some
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-14 01:48:30 UTC (rev 18365)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-14 03:08:08 UTC (rev 18366)
@@ -389,12 +389,9 @@
if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
call read_fault_files(localpath_name)
-
if (allocated(faults)) then
-
- call save_nodes_coords(nodes_coords,nnodes)
- call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
-
+ call save_nodes_coords(nodes_coords,nnodes)
+ call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
end if
end subroutine read_mesh_files
@@ -680,7 +677,7 @@
! endif
write(16,*) nnodes_loc
- call write_glob2loc_nodes_database(16, ipart, nnodes_loc, nodes_coords_virtual,&
+ call write_glob2loc_nodes_database(16, ipart, nnodes_loc, nodes_coords_open,&
glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, 2)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-14 01:48:30 UTC (rev 18365)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-14 03:08:08 UTC (rev 18366)
@@ -11,7 +11,7 @@
end type fault_type
type(fault_type), allocatable, save :: faults(:)
- double precision, dimension(:,:), allocatable, save :: nodes_coords_virtual
+ double precision, dimension(:,:), allocatable, save :: nodes_coords_open
integer, parameter :: long = SELECTED_INT_KIND(18)
@@ -20,7 +20,7 @@
! PGB: for a simple test is fine .For SCEC lower values.
public :: read_fault_files, fault_collect_elements, close_faults, write_fault_database, &
- save_nodes_coords, nodes_coords_virtual, faults
+ save_nodes_coords, nodes_coords_open, faults
CONTAINS
!==========================================================================================
@@ -114,8 +114,8 @@
integer, intent(in) :: nnodes
double precision, dimension(3,nnodes), intent(in) :: nodes_coords
- allocate(nodes_coords_virtual(3,nnodes))
- nodes_coords_virtual = nodes_coords
+ allocate(nodes_coords_open(3,nnodes))
+ nodes_coords_open = nodes_coords
end subroutine save_nodes_coords
@@ -179,15 +179,10 @@
dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
dist = sqrt(dist)
- !pgb : dist has to be > 0 in case we want to avoid re sticking nodes with
- !pgb : the same positions already , like nodes at fault edeges.
- !pgb : Or nodes that have been sticked already.
-
!jpa: Closing nodes that are already closed is not a problem
!jpa: I process them again to leave the loop as early as possible
!jpa: and to test if a facing node was found (See below).
- !if ((0.0d0 < dist) .and. (dist <= FAULT_GAP_TOLERANCE)) then
if (dist <= FAULT_GAP_TOLERANCE) then
xyz = (xyz_1 + xyz_2)*0.5d0
nodes_coords(:,iglob2) = xyz
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2011-05-14 01:48:30 UTC (rev 18365)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2011-05-14 03:08:08 UTC (rev 18366)
@@ -55,7 +55,7 @@
use create_regions_mesh_ext_par
use fault_object, only: fault_read_input,fault_setup, &
fault_save_arrays,fault_save_arrays_test,&
- nnodes_coords_open,nodes_coords_open,fault_db
+ nnodes_coords_open,nodes_coords_open,ANY_FAULT_IN_THIS_PROC
implicit none
@@ -169,66 +169,61 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
nspec2D_bottom,nspec2D_top,ANISOTROPY)
+ ! if faults exist this reads nodes_coords_open
+ call fault_read_input(prname,NDIM)
-! fills location and weights for Gauss-Lobatto-Legendre points, shape and derivations,
-! returns jacobianstore,xixstore,...gammazstore
-! and GLL-point locations in xstore,ystore,zstore
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...setting up jacobian '
endif
-! creates ibool index array for projection from local to global points
-
- call fault_read_input(prname,NDIM)
-
- if (allocated(fault_db)) then
-! Calculating jacobians with fault open to compute *store nedded for ibool.
- call crm_ext_setup_jacobian(myrank, &
+ if (ANY_FAULT_IN_THIS_PROC) then
+ ! compute jacobians with fault open and *store needed for ibool.
+ call crm_ext_setup_jacobian(myrank, &
xstore,ystore,zstore,nspec, &
nodes_coords_open, nnodes_coords_open,&
elmnts_ext_mesh,nelmnts_ext_mesh)
- call sync_all()
-! creates ibool index array for projection from local to global points
-! with fault open.
- if( myrank == 0) then
- write(IMAIN,*) ' ...indexing global points'
- endif
- call crm_ext_setup_indexing(ibool, &
+ ! create ibool with faults open
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...indexing global points'
+ endif
+ call crm_ext_setup_indexing(ibool, &
xstore,ystore,zstore,nspec,nglob,npointot, &
nnodes_coords_open,nodes_coords_open,myrank)
-! crm_ext_setup_indexing : computes xstore , ystore , zstore.
- call sync_all()
- if( myrank == 0) then
- write(IMAIN,*) ' ...setting up jacobian '
- endif
-! Sticking *store and calculating correctly jacobians.
- call crm_ext_setup_jacobian(myrank, &
+ ! recalculate *store with faults closed
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...setting up jacobian '
+ endif
+ call crm_ext_setup_jacobian(myrank, &
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
- else
-! With no fault.
- call crm_ext_setup_jacobian(myrank, &
+ ! at this point (xyz)store_dummy are still open
+
+ ! with no fault
+ else
+ call crm_ext_setup_jacobian(myrank, &
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
- call sync_all()
-! creates ibool index array for projection from local to global points
- if( myrank == 0) then
- write(IMAIN,*) ' ...indexing global points'
- endif
- call crm_ext_setup_indexing(ibool, &
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...indexing global points'
+ endif
+ call crm_ext_setup_indexing(ibool, &
xstore,ystore,zstore,nspec,nglob,npointot, &
nnodes_ext_mesh,nodes_coords_ext_mesh,myrank)
- end if
+ end if
call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
xstore,ystore,zstore,nspec,nglob,myrank)
-
+ ! this closes (xyz)store_dummy
+
! sets up MPI interfaces between partitions
call sync_all()
if( myrank == 0) then
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-14 01:48:30 UTC (rev 18365)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-14 03:08:08 UTC (rev 18366)
@@ -52,6 +52,7 @@
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
double precision, allocatable, save :: nodes_coords_open(:,:)
integer, save :: nnodes_coords_open
+ logical, save :: ANY_FAULT_IN_THIS_PROC = .false.
! corners indices of reference cube faces
integer,dimension(3,4),parameter :: iface1_corner_ijk = &
@@ -72,7 +73,7 @@
iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, fault_db, &
- nodes_coords_open, nnodes_coords_open
+ nodes_coords_open, nnodes_coords_open, ANY_FAULT_IN_THIS_PROC
contains
@@ -119,6 +120,8 @@
if (nspec == 0) cycle
+ ANY_FAULT_IN_THIS_PROC = .true.
+
allocate(fault_db(iflt)%ispec1(nspec))
allocate(fault_db(iflt)%inodes1(4,nspec))
allocate(fault_db(iflt)%ispec2(nspec))
@@ -165,21 +168,14 @@
double precision, dimension(NDIM,nnodes_ext_mesh),intent(in) :: nodes_coords_ext_mesh
integer :: iflt
- logical :: fault_exists
- if (.not. allocated(fault_db)) return
+ if (.not. ANY_FAULT_IN_THIS_PROC ) return
- fault_exists = .false.
do iflt=1,size(fault_db)
- if (fault_db(iflt)%nspec>0) then
- fault_exists = .true.
- exit
- endif
- enddo
- if (.not.fault_exists) return
-
- do iflt=1,size(fault_db)
+ ! close the fault in (xyz)store_dummy
+ call close_fault(fault_db(iflt))
+
call setup_iface(fault_db(iflt),nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool)
! saving gll indices for each fault face, needed for ibulks
@@ -195,13 +191,8 @@
call setup_Kelvin_Voigt_eta(fault_db(iflt),nspec)
- ! close the fault (modifies *store_dummy)
- call close_fault(fault_db(iflt))
-
call save_fault_xyzcoord_ibulk(fault_db(iflt),ibool,nspec)
- ! close the fault (modifies *store_dummy) to compute jacobians and normals
- call close_fault(fault_db(iflt)) !,xstore,ystore,zstore,nspec)
call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob,myrank)
enddo
@@ -210,7 +201,7 @@
!=============================================================================================================
-! We only close *store_dummy
+! We only close *store_dummy. Closing *store is more difficult.
! Fortunately only *store_dummy is needed to compute jacobians and normals
subroutine close_fault(fdb) !,xstore,ystore,zstore,nspec)
@@ -234,8 +225,10 @@
enddo
! Closing *store is not easy ...
-! the following would work only if ispec1 and ispec2 were always facing each other
-! and if igll corresponded to the same node in both element faces
+! The following would work only if ispec1 and ispec2 were always facing each other
+! and if igll corresponded to the same node in both element faces.
+! That would require reordering isepcs, e.g. sorting them by coordinate of the middle point of the face,
+! and then, within each face, searching the nearest node (by coordinate)
!
! do e=1,fdb%nspec
! e1 = fdb%ispec1(e)
@@ -596,16 +589,14 @@
integer :: nbfaults,iflt,ier
character(len=256) :: filename
+ if (.not.allocated(fault_db)) return
+
! saves mesh file proc***_fault_db.txt
filename = prname(1:len_trim(prname))//'fault_db.txt'
open(unit=IOUT,file=trim(filename),status='unknown',action='write',iostat=ier)
if( ier /= 0 ) stop 'error opening database proc######_external_mesh.bin'
- if (allocated(fault_db)) then
- nbfaults = size(fault_db)
- else
- nbfaults = 0
- endif
+ nbfaults = size(fault_db)
write(IOUT,*) 'NBFAULTS = ',nbfaults
do iflt=1,nbfaults
write(IOUT,*) 'BEGIN FAULT # ',iflt
@@ -671,6 +662,7 @@
integer :: nbfaults,iflt,ier
character(len=256) :: filename
+ if (.not.allocated(fault_db)) return
! saves mesh file proc***_Kelvin_voigt_eta.bin
if (allocated(Kelvin_Voigt_eta)) then
@@ -693,11 +685,7 @@
stop
endif
- if (allocated(fault_db)) then
- nbfaults = size(fault_db)
- else
- nbfaults = 0
- endif
+ nbfaults = size(fault_db)
write(IOUT) nbfaults
do iflt=1,nbfaults
call save_one_fault_bin(fault_db(iflt),IOUT)
More information about the CIG-COMMITS
mailing list