[cig-commits] r18362 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Fri May 13 15:18:09 PDT 2011
Author: percygalvez
Date: 2011-05-13 15:18:09 -0700 (Fri, 13 May 2011)
New Revision: 18362
Modified:
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:
updated
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-13 22:17:08 UTC (rev 18361)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2011-05-13 22:18:09 UTC (rev 18362)
@@ -53,8 +53,9 @@
! create the different regions of the mesh
use create_regions_mesh_ext_par
- use fault_object, only: fault_read_input, fault_setup, &
- fault_save_arrays , fault_save_arrays_test
+ use fault_object, only: fault_read_input,fault_setup, &
+ fault_save_arrays,fault_save_arrays_test,&
+ nnodes_coords_open,nodes_coords_open,fault_db
implicit none
@@ -177,25 +178,54 @@
write(IMAIN,*) ' ...setting up jacobian '
endif
-!jpa: jacobians should be computed after closgin the fault (after fault_setup)
-!jpa: can we tolerate this small error?
- call crm_ext_setup_jacobian(myrank, &
+! 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, &
+ 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, &
+ 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, &
+ 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, &
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
-! creates ibool index array for projection from local to global points
- call sync_all()
- if( myrank == 0) then
- write(IMAIN,*) ' ...indexing global points'
- endif
-
-! crm_ext_setup_indexing : computes xstore , ystore , zstore.
- call crm_ext_setup_indexing(ibool, &
+ 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, &
xstore,ystore,zstore,nspec,nglob,npointot, &
nnodes_ext_mesh,nodes_coords_ext_mesh,myrank)
+ end if
- call fault_read_input(prname)
call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
xstore,ystore,zstore,nspec,nglob,myrank)
@@ -608,7 +638,6 @@
!
!-------------------------------------------------------------------------------------------------
-!
subroutine crm_ext_setup_indexing(ibool, &
xstore,ystore,zstore,nspec,nglob,npointot, &
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-13 22:17:08 UTC (rev 18361)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-13 22:18:09 UTC (rev 18362)
@@ -50,7 +50,9 @@
type(fault_db_type), allocatable, save :: fault_db(:)
! fault_db(i) is the database of the i-th fault in the mesh
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
-
+ double precision, allocatable, save :: nodes_coords_open(:,:)
+ integer, save :: nnodes_coords_open
+
! 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
@@ -69,19 +71,21 @@
iface3_corner_ijk,iface4_corner_ijk, &
iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays
+ public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, fault_db, &
+ nodes_coords_open, nnodes_coords_open
contains
!=================================================================================================================
-subroutine fault_read_input(prname)
+subroutine fault_read_input(prname,NDIM)
character(len=256), intent(in) :: prname
-
- integer :: nb,i,iflt,ier,nspec
+ integer, intent(in) :: NDIM
+
+ integer :: nb,i,iflt,ier,nspec,dummy_node
integer, parameter :: IIN = 100
real(kind=CUSTOM_REAL) :: eta
-
+
! read fault input file
nb = 0
open(unit=IIN,file='DATA/FAULT/Par_file_faults.in',status='old',action='read',iostat=ier)
@@ -122,13 +126,11 @@
do i=1,nspec
read(IIN,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%inodes1(:,i)
- write(6,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%inodes1(:,i) !TEST
enddo
do i=1,nspec
read(IIN,*) fault_db(iflt)%ispec2(i), fault_db(iflt)%inodes2(:,i)
- write(6,*) fault_db(iflt)%ispec2(i), fault_db(iflt)%inodes2(:,i) !TEST
enddo
-
+
! loading ispec1 ispec2 iface1 iface2 of fault elements.
! allocate(fault_db(iflt)%iface1(nspec))
! allocate(fault_db(iflt)%iface2(nspec))
@@ -139,6 +141,13 @@
enddo
+ ! read nodes_coords with fault crack open.
+ read(IIN,*) nnodes_coords_open
+ allocate(nodes_coords_open(NDIM,nnodes_coords_open))
+ do i = 1, nnodes_coords_open
+ read(IIN,*) dummy_node, nodes_coords_open(:,i)
+ enddo
+
close(IIN)
end subroutine fault_read_input
@@ -185,7 +194,12 @@
call setup_ibulks(fault_db(iflt),ibool,nspec)
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)
@@ -194,6 +208,7 @@
end subroutine fault_setup
+
!=============================================================================================================
! We only close *store_dummy
! Fortunately only *store_dummy is needed to compute jacobians and normals
@@ -407,8 +422,6 @@
subroutine setup_ibulks(fdb,ibool,nspec)
- use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
-
type(fault_db_type), intent(inout) :: fdb
integer, intent(in) :: nspec, ibool(NGLLX,NGLLY,NGLLZ,nspec)
@@ -416,13 +429,6 @@
allocate( fdb%ibulk1(fdb%nglob) )
allocate( fdb%ibulk2(fdb%nglob) )
- allocate( fdb%xcoordbulk1(fdb%nglob) )
- allocate( fdb%ycoordbulk1(fdb%nglob) )
- allocate( fdb%zcoordbulk1(fdb%nglob) )
- allocate( fdb%xcoordbulk2(fdb%nglob) )
- allocate( fdb%ycoordbulk2(fdb%nglob) )
- allocate( fdb%zcoordbulk2(fdb%nglob) )
-
do e=1, fdb%nspec
do k=1, NGLLSQUARE
@@ -432,27 +438,78 @@
ke=fdb%ijk1(3,k,e)
K1= fdb%ibool1(k,e)
fdb%ibulk1(K1)=ibool(ie,je,ke,fdb%ispec1(e))
- fdb%xcoordbulk1(K1) = xstore_dummy(fdb%ibulk1(K1))
- fdb%ycoordbulk1(K1) = ystore_dummy(fdb%ibulk1(K1))
- fdb%zcoordbulk1(K1) = zstore_dummy(fdb%ibulk1(K1))
ie=fdb%ijk2(1,k,e)
je=fdb%ijk2(2,k,e)
ke=fdb%ijk2(3,k,e)
K2= fdb%ibool2(k,e)
fdb%ibulk2(K2)=ibool(ie,je,ke,fdb%ispec2(e))
- fdb%xcoordbulk2(K2) = xstore_dummy(fdb%ibulk2(K2))
- fdb%ycoordbulk2(K2) = ystore_dummy(fdb%ibulk2(K2))
- fdb%zcoordbulk2(K2) = zstore_dummy(fdb%ibulk2(K2))
-
+
enddo
enddo
end subroutine setup_ibulks
+!=============================================================================================================
+! We only close *store_dummy. *store is close already in create_regions_mesh.f90
+! Fortunately only *store_dummy is needed to compute jacobians and normals
+subroutine close_fault(fdb)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
+
+ type(fault_db_type), intent(inout) :: fdb
+
+ integer :: i,K1,K2
+
+ do i=1,fdb%nglob
+ K1 = fdb%ibulk1(i)
+ K2 = fdb%ibulk2(i)
+ xstore_dummy(K1) = 0.5d0*( xstore_dummy(K1) + xstore_dummy(K2) )
+ xstore_dummy(K2) = xstore_dummy(K1)
+ ystore_dummy(K1) = 0.5d0*( ystore_dummy(K1) + ystore_dummy(K2) )
+ ystore_dummy(K2) = ystore_dummy(K1)
+ zstore_dummy(K1) = 0.5d0*( zstore_dummy(K1) + zstore_dummy(K2) )
+ zstore_dummy(K2) = zstore_dummy(K1)
+ enddo
+
+end subroutine close_fault
+
!=================================================================================
+ subroutine save_fault_xyzcoord_ibulk(fdb,ibool,nspec)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
+
+ type(fault_db_type), intent(inout) :: fdb
+ integer, intent(in) :: nspec, ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer :: K1, K2, i
+
+ allocate( fdb%xcoordbulk1(fdb%nglob) )
+ allocate( fdb%ycoordbulk1(fdb%nglob) )
+ allocate( fdb%zcoordbulk1(fdb%nglob) )
+ allocate( fdb%xcoordbulk2(fdb%nglob) )
+ allocate( fdb%ycoordbulk2(fdb%nglob) )
+ allocate( fdb%zcoordbulk2(fdb%nglob) )
+
+ do i=1, fdb%nglob
+ K1 =fdb%ibulk1(i)
+ K2 =fdb%ibulk2(i)
+ fdb%xcoordbulk1(i) = xstore_dummy(K1)
+ fdb%ycoordbulk1(i) = ystore_dummy(K1)
+ fdb%zcoordbulk1(i) = zstore_dummy(K1)
+
+ fdb%xcoordbulk2(i) = xstore_dummy(K2)
+ fdb%ycoordbulk2(i) = ystore_dummy(K2)
+ fdb%zcoordbulk2(i) = zstore_dummy(K2)
+ enddo
+
+ end subroutine save_fault_xyzcoord_ibulk
+
+
+!=================================================================================
+
subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob,myrank)
use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy, &
More information about the CIG-COMMITS
mailing list