[cig-commits] r21109 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh generate_databases specfem3D
surendra at geodynamics.org
surendra at geodynamics.org
Thu Dec 6 20:28:14 PST 2012
Author: surendra
Date: 2012-12-06 20:28:13 -0800 (Thu, 06 Dec 2012)
New Revision: 21109
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/fault_generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_global.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90
Log:
Fixed some compile time errors
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/Makefile.in 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/Makefile.in 2012-12-07 04:28:13 UTC (rev 21109)
@@ -131,12 +131,12 @@
$O/part_decompose_mesh.o: part_decompose_mesh.f90
${FCCOMPILE_CHECK} -c -o $O/part_decompose_mesh.o part_decompose_mesh.f90 $(SCOTCH_INC)
-$O/decompose_mesh.o: decompose_mesh.F90 part_decompose_mesh.f90 $O/part_decompose_mesh.o
+$O/fault_scotch.o: fault_scotch.f90
+ ${FCCOMPILE_CHECK} -c -o $O/fault_scotch.o fault_scotch.f90 $(SCOTCH_INC)
+
+$O/decompose_mesh.o: decompose_mesh.F90 part_decompose_mesh.f90 $O/part_decompose_mesh.o $O/fault_scotch.o
${FCCOMPILE_CHECK} -c -o $O/decompose_mesh.o decompose_mesh.F90 $(SCOTCH_INC)
-$O/fault_scotch.o: fault_scotch.f90 part_decompose_mesh.f90 $O/part_decompose_mesh.o
- ${FCCOMPILE_CHECK} -c -o $O/fault_scotch.o fault_scotch.f90 $(SCOTCH_INC)
-
$O/program_decompose_mesh.o: program_decompose_mesh.f90 $O/part_decompose_mesh.o $O/decompose_mesh.o
${FCCOMPILE_CHECK} -c -o $O/program_decompose_mesh.o program_decompose_mesh.f90 $(SCOTCH_INC)
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -38,6 +38,7 @@
module decompose_mesh
use part_decompose_mesh
+ use fault_scotch
implicit none
@@ -866,7 +867,7 @@
deallocate(num_material)
! re-partitioning transfers two coupled elements on fault side 1 and side 2 to the same partition
- if (ANY_FAULT) call fault_repartition (nspec, nnodes, elmnts, nsize, nparts, part, esize, nodes_coords)
+ if (ANY_FAULT) call fault_repartition (nspec, nnodes, elmnts, nsize, nparts, part, NGNOD, nodes_coords)
! re-partitioning puts moho-surface coupled elements into same partition
@@ -1016,6 +1017,30 @@
close(IIN_database)
+
+ ! write fault database
+ if (ANY_FAULT) then
+ write(prname, "(i6.6,'_Database_fault')") ipart
+ open(unit=16,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
+ status='replace', action='write', form='unformatted', iostat = ier)
+ if( ier /= 0 ) then
+ print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
+ print*
+ print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
+ stop
+ endif
+ call write_fault_database(16, ipart, nspec, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, part)
+ !write(16,*) nnodes_loc
+ write(16) nnodes_loc
+ call write_glob2loc_nodes_database(16, ipart, nnodes_loc, nodes_coords_open,&
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 2)
+ close(16)
+ endif
+
+
enddo
print*, 'partitions: '
print*, ' num = ',nparts
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/fault_scotch.f90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -15,7 +15,7 @@
double precision, dimension(:,:), allocatable, save :: nodes_coords_open
logical, save :: ANY_FAULT = .false.
- logical, parameter :: PARALLEL_FAULT = .false.
+ logical, parameter :: PARALLEL_FAULT = .true.
! NOTE: PARALLEL_FAULT has to be the same
! in fault_solver_common.f90, fault_generate_databases.f90 and fault_scotch.f90
@@ -347,7 +347,7 @@
subroutine fault_repartition (nelmnts, nnodes, elmnts, nsize, nproc, part, esize, nodes_coords)
- integer(long), intent(in) :: nelmnts,nsize
+ integer, intent(in) :: nelmnts,nsize
integer, intent(in) :: nnodes, nproc, esize
integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
integer, dimension(0:nelmnts-1), intent(inout) :: part
@@ -519,7 +519,8 @@
print *, 'Fatal error: Number of fault elements do not coincide. Abort.'
stop
end if
- write(IIN_database,*) nspec_fault_1
+ !write(IIN_database,*) nspec_fault_1
+ write(IIN_database) nspec_fault_1
! if no fault element in this partition, move to next fault
if (nspec_fault_1==0) cycle
@@ -536,7 +537,8 @@
end if
end do
end do
- write(IIN_database,*) glob2loc_elmnts(e-1)+1, loc_nodes
+ !write(IIN_database,*) glob2loc_elmnts(e-1)+1, loc_nodes
+ write(IIN_database) glob2loc_elmnts(e-1)+1, loc_nodes
end if
enddo
@@ -552,7 +554,8 @@
end if
end do
end do
- write(IIN_database,*) glob2loc_elmnts(e-1)+1, loc_nodes
+ !write(IIN_database,*) glob2loc_elmnts(e-1)+1, loc_nodes
+ write(IIN_database) glob2loc_elmnts(e-1)+1, loc_nodes
end if
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/Makefile.in 2012-12-07 04:28:13 UTC (rev 21109)
@@ -70,6 +70,7 @@
$O/tomography_par.o \
$O/assemble_MPI_scalar.shared.o \
$O/calc_jacobian.o \
+ $O/fault_generate_databases.o \
$O/check_mesh_resolution.shared.o \
$O/create_name_database.shared.o \
$O/create_mass_matrices.o \
@@ -78,7 +79,6 @@
$O/define_derivation_matrices.shared.o \
$O/detect_surface.shared.o \
$O/exit_mpi.shared.o \
- $O/fault_generate_databases.o \
$O/finalize_databases.o \
$O/generate_databases.o \
$O/get_absorbing_boundary.o \
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -51,6 +51,11 @@
nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho
use create_regions_mesh_ext_par
+ use fault_generate_databases, only: fault_read_input,fault_setup, &
+ fault_save_arrays,fault_save_arrays_test,&
+ nnodes_coords_open,nodes_coords_open,ANY_FAULT_IN_THIS_PROC,&
+ ANY_FAULT, PARALLEL_FAULT
+
implicit none
! local parameters
@@ -68,6 +73,9 @@
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,myrank)
+
! fills location and weights for Gauss-Lobatto-Legendre points, shape and derivations,
! returns jacobianstore,xixstore,...gammazstore
! and GLL-point locations in xstore,ystore,zstore
@@ -75,20 +83,50 @@
if( myrank == 0) then
write(IMAIN,*) ' ...setting up jacobian '
endif
- 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)
+ 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)
+ endif
+
! creates ibool index array for projection from local to global points
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)
+ if (ANY_FAULT_IN_THIS_PROC) then
+ call crm_ext_setup_indexing(ibool, &
+ xstore,ystore,zstore,nspec,nglob,npointot, &
+ nnodes_coords_open,nodes_coords_open,myrank)
+ else ! with no fault
+ call crm_ext_setup_indexing(ibool, &
+ xstore,ystore,zstore,nspec,nglob,npointot, &
+ nnodes_ext_mesh,nodes_coords_ext_mesh,myrank)
+ end if
+ if (ANY_FAULT) then
+ ! recalculate *store with faults closed
+ call sync_all()
+ if (myrank == 0) write(IMAIN,*) ' ... resetting up jacobian in fault domains'
+ if (ANY_FAULT_IN_THIS_PROC) call crm_ext_setup_jacobian(myrank, &
+ xstore,ystore,zstore,nspec, &
+ nodes_coords_ext_mesh,nnodes_ext_mesh,&
+ elmnts_ext_mesh,nelmnts_ext_mesh)
+ ! at this point (xyz)store_dummy are still open
+ if (.NOT. PARALLEL_FAULT) call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
+ xstore,ystore,zstore,nspec,nglob,myrank)
+ ! this closes (xyz)store_dummy
+ endif
+
+
! sets up MPI interfaces between partitions
call sync_all()
if( myrank == 0) then
@@ -102,6 +140,17 @@
num_interfaces_ext_mesh,max_interface_size_ext_mesh,&
my_neighbours_ext_mesh)
+
+ !SURENDRA (setting up parallel fault)
+ if (PARALLEL_FAULT .AND. ANY_FAULT) then
+ call sync_all()
+ !at this point (xyz)store_dummy are still open
+ call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
+ xstore,ystore,zstore,nspec,nglob,myrank)
+ ! this closes (xyz)store_dummy
+ endif
+
+
! sets up absorbing/free surface boundaries
call sync_all()
if( myrank == 0) then
@@ -183,6 +232,9 @@
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
SAVE_MESH_FILES,ANISOTROPY)
+! call fault_save_arrays_test(prname) ! for debugging
+ call fault_save_arrays(prname)
+
! saves moho surface
if( SAVE_MOHO_MESH ) then
call crm_save_moho()
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/fault_generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/fault_generate_databases.f90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/fault_generate_databases.f90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -9,7 +9,8 @@
module fault_generate_databases
- use create_regions_mesh_ext_par, only: NGLLX,NGLLY,NGLLZ,NGLLSQUARE,NGNOD2D,NDIM,CUSTOM_REAL,IMAIN
+ use create_regions_mesh_ext_par, only: NGLLX,NGLLY,NGLLZ,NGLLSQUARE,NDIM,CUSTOM_REAL,IMAIN
+ use generate_databases_par, only : NGNOD2D
implicit none
private
@@ -72,7 +73,7 @@
! read fault input file
nb = 0
- open(unit=IIN_PAR,file='DATA/Par_file_faults',status='old',action='read',iostat=ier)
+ open(unit=IIN_PAR,file='../DATA/Par_file_faults',status='old',action='read',iostat=ier)
if (ier==0) then
read(IIN_PAR,*) nb
if (myrank==0) write(IMAIN,*) ' ... reading ', nb,' faults from file DATA/Par_file_faults'
@@ -92,16 +93,16 @@
! read fault database file
open(unit=IIN_PAR,file=prname(1:len_trim(prname))//'Database_fault', &
- status='old',action='read',form='formatted',iostat=ier)
+ status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
- write(IIN_PAR,*) 'error opening file: ',prname(1:len_trim(prname))//'Database_fault'
- write(IIN_PAR,*) 'make sure file exists'
+ write(IMAIN,*) 'error opening file: ',prname(1:len_trim(prname))//'Database_fault'
+ write(IMAIN,*) 'make sure file exists'
stop
endif
do iflt=1,size(fault_db)
- read(IIN_PAR,*) nspec
+ read(IIN_PAR) nspec
fault_db(iflt)%nspec = nspec
if (nspec == 0) cycle
@@ -114,10 +115,10 @@
allocate(fault_db(iflt)%inodes2(4,nspec))
do i=1,nspec
- read(IIN_PAR,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%inodes1(:,i)
+ read(IIN_PAR) fault_db(iflt)%ispec1(i), fault_db(iflt)%inodes1(:,i)
enddo
do i=1,nspec
- read(IIN_PAR,*) fault_db(iflt)%ispec2(i), fault_db(iflt)%inodes2(:,i)
+ read(IIN_PAR) fault_db(iflt)%ispec2(i), fault_db(iflt)%inodes2(:,i)
enddo
! loading ispec1 ispec2 iface1 iface2 of fault elements.
@@ -131,10 +132,10 @@
enddo
! read nodes coordinates of the original version of the mesh, in which faults are open
- read(IIN_PAR,*) nnodes_coords_open
+ read(IIN_PAR) nnodes_coords_open
allocate(nodes_coords_open(NDIM,nnodes_coords_open))
do i = 1, nnodes_coords_open
- read(IIN_PAR,*) dummy_node, nodes_coords_open(:,i)
+ read(IIN_PAR) dummy_node, nodes_coords_open(:,i)
enddo
close(IIN_PAR)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_global.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_global.f90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_global.f90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -55,14 +55,12 @@
! small value for double precision and to avoid sensitivity to roundoff
double precision SMALLVALTOL
-! number of points per spectral element
-! integer, parameter :: NGLLCUBE = NGLLX * NGLLY * NGLLZ
!jpampuero To allow usage of this routine for volume and surface meshes:
!jpampuero For volumes NGLLCUBE = NGLLX * NGLLY * NGLLZ
!jpampuero For surfaces NGLLCUBE = NGLLX * NGLLY
- integer :: NGLLCUBE
+ integer :: NGLLCUBE_local
- NGLLCUBE=npointot/nspec
+ NGLLCUBE_local=npointot/nspec
! for vectorization of loops
! integer, parameter :: NGLLCUBE_NDIM = NGLLCUBE * NDIM
@@ -78,8 +76,8 @@
! establish initial pointers
do ispec=1,nspec
- ieoff=NGLLCUBE*(ispec-1)
- do ilocnum=1,NGLLCUBE
+ ieoff=NGLLCUBE_local*(ispec-1)
+ do ilocnum=1,NGLLCUBE_local
loc(ilocnum+ieoff)=ilocnum+ieoff
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2012-12-07 04:28:13 UTC (rev 21109)
@@ -114,9 +114,6 @@
$O/define_derivation_matrices.shared.o \
$O/detect_surface.shared.o \
$O/exit_mpi.shared.o \
- $O/fault_solver_common.o \
- $O/fault_solver_dynamic.o \
- $O/fault_solver_kinematic.o \
$O/force_ftz.cc.o \
$O/get_attenuation_model.shared.o \
$O/get_cmt.shared.o \
@@ -178,6 +175,9 @@
$O/specfem3D_par.o \
$O/PML_init.o \
$O/assemble_MPI_vector.o \
+ $O/fault_solver_common.o \
+ $O/fault_solver_dynamic.o \
+ $O/fault_solver_kinematic.o \
$O/compute_add_sources_acoustic.o \
$O/compute_add_sources_elastic.o \
$O/compute_add_sources_poroelastic.o \
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic_noDev.f90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -175,7 +175,9 @@
integer ispec,iglob,ispec_p,num_elements
integer i,j,k,l
+ real(kind=CUSTOM_REAL) :: eta
+
if( iphase == 1 ) then
num_elements = nspec_outer_elastic
else
@@ -239,23 +241,23 @@
do l=1,NGLLX
hp1 = hprime_xx(i,l)
iglob = ibool(l,j,k,ispec)
- tempx1l = tempx1l + dloc(1,iglob)*hp1
- tempy1l = tempy1l + dloc(2,iglob)*hp1
- tempz1l = tempz1l + dloc(3,iglob)*hp1
+ tempx1l = tempx1l + dloc(1,l,j,k)*hp1
+ tempy1l = tempy1l + dloc(2,l,j,k)*hp1
+ tempz1l = tempz1l + dloc(3,l,j,k)*hp1
!!! can merge these loops because NGLLX = NGLLY = NGLLZ
hp2 = hprime_yy(j,l)
iglob = ibool(i,l,k,ispec)
- tempx2l = tempx2l + dloc(1,iglob)*hp2
- tempy2l = tempy2l + dloc(2,iglob)*hp2
- tempz2l = tempz2l + dloc(3,iglob)*hp2
+ tempx2l = tempx2l + dloc(1,i,l,k)*hp2
+ tempy2l = tempy2l + dloc(2,i,l,k)*hp2
+ tempz2l = tempz2l + dloc(3,i,l,k)*hp2
!!! can merge these loops because NGLLX = NGLLY = NGLLZ
hp3 = hprime_zz(k,l)
iglob = ibool(i,j,l,ispec)
- tempx3l = tempx3l + dloc(1,iglob)*hp3
- tempy3l = tempy3l + dloc(2,iglob)*hp3
- tempz3l = tempz3l + dloc(3,iglob)*hp3
+ tempx3l = tempx3l + dloc(1,i,j,l)*hp3
+ tempy3l = tempy3l + dloc(2,i,j,l)*hp3
+ tempz3l = tempz3l + dloc(3,i,j,l)*hp3
enddo
if( ATTENUATION .and. COMPUTE_AND_STORE_STRAIN ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90 2012-12-06 09:07:32 UTC (rev 21108)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90 2012-12-07 04:28:13 UTC (rev 21109)
@@ -320,9 +320,12 @@
real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
real(kind=CUSTOM_REAL) :: r1(size(a))
integer :: i
+ real(kind=CUSTOM_REAL) :: SMALLVAL
NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+ SMALLVAL = 1.e-10_CUSTOM_REAL
+
if (n==0) return
do i=1,n
More information about the CIG-COMMITS
mailing list