[cig-commits] r19245 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh_SCOTCH decompose_mesh_SCOTCH/scotch_5.1.11/src/libscotch decompose_mesh_SCOTCH/scotch_5.1.11/src/scotch generate_databases shared specfem3D
cmorency at geodynamics.org
cmorency at geodynamics.org
Mon Nov 28 20:02:33 PST 2011
Author: cmorency
Date: 2011-11-28 20:02:32 -0800 (Mon, 28 Nov 2011)
New Revision: 19245
Added:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/libscotch/Makefile
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/scotch/Makefile
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
Log:
Merging of 3D Biot into SPECFEM3D with CUBIT support: corrected some bu with acoustic-poroelastic coupling. Added elastic-poroelastic coupling.
Note: coupling surfaces search for acoustic-poroelastic needs to identify poroelastic elements (properties needed for coupling). coupling surfaces search for poroelastic-elastic needs to identify both poroelastic and elastic elements (properties are needed and tractions are been calculated).
For now coupled elastic-poroelastic is implemented without Deville optimizations - to come soon.
Still in progress: poroelastic adjoint.
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -699,7 +699,7 @@
if( ier /= 0 ) stop 'error allocating array num_material'
num_material(:) = mat(1,:)
- ! in case of acoustic/elastic simulation, weights elements accordingly
+ ! in case of acoustic/elastic/poro simulation, weights elements accordingly
call acoustic_elastic_poro_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
num_material,mat_prop,undef_mat_prop)
@@ -790,13 +790,14 @@
endif
- ! re-partitioning puts acoustic-elastic coupled elements into same partition
+ ! re-partitioning puts poroelastic-elastic coupled elements into same partition
! integer :: nfaces_coupled
! integer, dimension(:,:), pointer :: faces_coupled
- ! call acoustic_elastic_repartitioning (nspec, nnodes, elmnts, &
- ! count_def_mat, mat(1,:) , mat_prop, &
- ! sup_neighbour, nsize, &
- ! nparts, part, nfaces_coupled, faces_coupled)
+ call poroelastic_elastic_repartitioning (nspec, nnodes, elmnts, &
+ count_def_mat, mat(1,:) , mat_prop, &
+ sup_neighbour, nsize, &
+ nparts, part)
+ !nparts, part, nfaces_coupled, faces_coupled)
! re-partitioning puts moho-surface coupled elements into same partition
call moho_surface_repartitioning (nspec, nnodes, elmnts, &
@@ -812,7 +813,7 @@
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, nparts)
! mpi interfaces
- ! acoustic/elastic boundaries will be split into different MPI partitions
+ ! acoustic/elastic/poroelastic boundaries will be split into different MPI partitions
call build_interfaces(nspec, sup_neighbour, part, elmnts, &
xadj, adjncy, tab_interfaces, &
tab_size_interfaces, ninterfaces, &
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -1372,13 +1372,14 @@
!--------------------------------------------------
- ! Repartitioning : two coupled acoustic/elastic elements are transfered to the same partition
+ ! Repartitioning : two coupled poroelastic/elastic elements are transfered to the same partition
!--------------------------------------------------
- subroutine acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, &
+ subroutine poroelastic_elastic_repartitioning (nelmnts, nnodes, elmnts, &
nb_materials, num_material, mat_prop, &
sup_neighbour, nsize, &
- nproc, part, nfaces_coupled, faces_coupled)
+ nproc, part)
+ !nproc, part, nfaces_coupled, faces_coupled)
implicit none
@@ -1393,11 +1394,12 @@
integer, dimension(0:nelmnts-1) :: part
integer, dimension(0:esize*nelmnts-1) :: elmnts
- integer, intent(out) :: nfaces_coupled
+ integer :: nfaces_coupled
+ !integer, intent(out) :: nfaces_coupled
integer, dimension(:,:), pointer :: faces_coupled
- logical, dimension(nb_materials) :: is_acoustic, is_elastic
+ logical, dimension(nb_materials) :: is_poroelastic, is_elastic
! neighbors
integer, dimension(:), allocatable :: xadj
@@ -1410,12 +1412,12 @@
integer :: el, el_adj
logical :: is_repartitioned
- ! sets acoustic/elastic flags for materials
- is_acoustic(:) = .false.
+ ! sets poroelastic/elastic flags for materials
+ is_poroelastic(:) = .false.
is_elastic(:) = .false.
do i = 1, nb_materials
- if (mat_prop(6,i) == 1 ) then
- is_acoustic(i) = .true.
+ if (mat_prop(6,i) == 3 ) then
+ is_poroelastic(i) = .true.
endif
if (mat_prop(6,i) == 2 ) then
is_elastic(i) = .true.
@@ -1439,7 +1441,7 @@
! counts coupled elements
nfaces_coupled = 0
do el = 0, nelmnts-1
- if ( is_acoustic(num_material(el+1)) ) then
+ if ( is_poroelastic(num_material(el+1)) ) then
do el_adj = xadj(el), xadj(el+1) - 1
if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
nfaces_coupled = nfaces_coupled + 1
@@ -1455,7 +1457,7 @@
! stores elements indices
nfaces_coupled = 0
do el = 0, nelmnts-1
- if ( is_acoustic(num_material(el+1)) ) then
+ if ( is_poroelastic(num_material(el+1)) ) then
do el_adj = xadj(el), xadj(el+1) - 1
if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
nfaces_coupled = nfaces_coupled + 1
@@ -1484,7 +1486,7 @@
endif
enddo
- end subroutine acoustic_elastic_repartitioning
+ end subroutine poroelastic_elastic_repartitioning
!--------------------------------------------------
! Repartitioning : two coupled moho surface elements are transfered to the same partition
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/libscotch/Makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/libscotch/Makefile 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/libscotch/Makefile 2011-11-29 04:02:32 UTC (rev 19245)
@@ -97,7 +97,20 @@
$(libdir)/libscotcherrexit$(LIB) \
$(includedir)/scotch.h $(includedir)/scotchf.h
+$(libdir)/libscotcherr$(LIB): libscotcherr$(LIB)
+ -$(CP) libscotcherr$(LIB) $(libdir)
+$(libdir)/libscotcherrexit$(LIB): libscotcherrexit$(LIB)
+ -$(CP) libscotcherrexit$(LIB) $(libdir)
+$(includedir)/scotch.h: scotch.h
+ -$(CP) scotch.h $(includedir)
+$(includedir)/scotchf.h: scotchf.h
+ -$(CP) scotchf.h $(includedir)
+install : $(libdir)/libscotch$(LIB) $(libdir)/libscotcherr$(LIB) \
+ $(libdir)/libscotcherrexit$(LIB) \
+ $(includedir)/scotch.h $(includedir)/scotchf.h
+
+
ptinstall :
-$(CP) scotch.h $(includedir)/ptscotch.h
-$(CP) scotchf.h $(includedir)/ptscotchf.h
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/scotch/Makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/scotch/Makefile 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/scotch_5.1.11/src/scotch/Makefile 2011-11-29 04:02:32 UTC (rev 19245)
@@ -241,7 +241,78 @@
$(bindir)/mmk_m2$(EXE) $(bindir)/mmk_m3$(EXE) \
$(bindir)/mord$(EXE) $(bindir)/mtst$(EXE)
+ -$(CP) amk_ccc$(EXE) $(bindir)/amk_ccc$(EXE)
+$(bindir)/amk_fft2$(EXE): amk_fft2$(EXE)
+ -$(CP) amk_fft2$(EXE) $(bindir)/amk_fft2$(EXE)
+$(bindir)/amk_grf$(EXE): amk_grf$(EXE)
+ -$(CP) amk_grf$(EXE) $(bindir)/amk_grf$(EXE)
+$(bindir)/amk_hy$(EXE): amk_hy$(EXE)
+ -$(CP) amk_hy$(EXE) $(bindir)/amk_hy$(EXE)
+$(bindir)/amk_m2$(EXE): amk_m2$(EXE)
+ -$(CP) amk_m2$(EXE) $(bindir)/amk_m2$(EXE)
+$(bindir)/amk_p2$(EXE): amk_p2$(EXE)
+ -$(CP) amk_p2$(EXE) $(bindir)/amk_p2$(EXE)
+$(bindir)/atst$(EXE): atst$(EXE)
+ -$(CP) atst$(EXE) $(bindir)/atst$(EXE)
+$(bindir)/gbase$(EXE): gbase$(EXE)
+ -$(CP) gbase$(EXE) $(bindir)/gbase$(EXE)
+$(bindir)/gcv$(EXE): gcv$(EXE)
+ -$(CP) gcv$(EXE) $(bindir)/gcv$(EXE)
+$(bindir)/gmap$(EXE): gmap$(EXE)
+ -$(CP) gmap$(EXE) $(bindir)/gmap$(EXE)
+$(bindir)/gmk_hy$(EXE): gmk_hy$(EXE)
+ -$(CP) gmk_hy$(EXE) $(bindir)/gmk_hy$(EXE)
+$(bindir)/gmk_m2$(EXE): gmk_m2$(EXE)
+ -$(CP) gmk_m2$(EXE) $(bindir)/gmk_m2$(EXE)
+$(bindir)/gmk_m3$(EXE): gmk_m3$(EXE)
+ -$(CP) gmk_m3$(EXE) $(bindir)/gmk_m3$(EXE)
+$(bindir)/gmk_msh$(EXE): gmk_msh$(EXE)
+ -$(CP) gmk_msh$(EXE) $(bindir)/gmk_msh$(EXE)
+$(bindir)/gmk_ub2$(EXE): gmk_ub2$(EXE)
+ -$(CP) gmk_ub2$(EXE) $(bindir)/gmk_ub2$(EXE)
+$(bindir)/gmtst$(EXE): gmtst$(EXE)
+ -$(CP) gmtst$(EXE) $(bindir)/gmtst$(EXE)
+$(bindir)/gord$(EXE): gord$(EXE)
+ -$(CP) gord$(EXE) $(bindir)/gord$(EXE)
+$(bindir)/gotst$(EXE): gotst$(EXE)
+ -$(CP) gotst$(EXE) $(bindir)/gotst$(EXE)
+$(bindir)/gout$(EXE): gout$(EXE)
+ -$(CP) gout$(EXE) $(bindir)/gout$(EXE)
+$(bindir)/gpart$(EXE): $(bindir)/gmap$(EXE)
+ -$(RM) $(bindir)/gpart$(EXE)
+ -$(LN) $(bindir)/gmap$(EXE) $(bindir)/gpart$(EXE)
+$(bindir)/gscat$(EXE): gscat$(EXE)
+ -$(CP) gscat$(EXE) $(bindir)/gscat$(EXE)
+$(bindir)/gtst$(EXE): gtst$(EXE)
+ -$(CP) gtst$(EXE) $(bindir)/gtst$(EXE)
+$(bindir)/mcv$(EXE): mcv$(EXE)
+ -$(CP) mcv$(EXE) $(bindir)/mcv$(EXE)
+$(bindir)/mmk_m2$(EXE): mmk_m2$(EXE)
+ -$(CP) mmk_m2$(EXE) $(bindir)/mmk_m2$(EXE)
+$(bindir)/mmk_m3$(EXE): mmk_m3$(EXE)
+ -$(CP) mmk_m3$(EXE) $(bindir)/mmk_m3$(EXE)
+$(bindir)/mord$(EXE): mord$(EXE)
+ -$(CP) mord$(EXE) $(bindir)/mord$(EXE)
+$(bindir)/mtst$(EXE): mtst$(EXE)
+ -$(CP) mtst$(EXE) $(bindir)/mtst$(EXE)
+
+install : $(bindir)/acpl$(EXE) $(bindir)/amk_ccc$(EXE) \
+ $(bindir)/amk_fft2$(EXE) $(bindir)/amk_grf$(EXE) \
+ $(bindir)/amk_hy$(EXE) $(bindir)/amk_m2$(EXE) \
+ $(bindir)/amk_p2$(EXE) $(bindir)/atst$(EXE) \
+ $(bindir)/gbase$(EXE) $(bindir)/gcv$(EXE) \
+ $(bindir)/gmap$(EXE) $(bindir)/gmk_hy$(EXE) \
+ $(bindir)/gmk_m2$(EXE) $(bindir)/gmk_m3$(EXE) \
+ $(bindir)/gmk_msh$(EXE) $(bindir)/gmk_ub2$(EXE) \
+ $(bindir)/gmtst$(EXE) $(bindir)/gord$(EXE) \
+ $(bindir)/gotst$(EXE) $(bindir)/gout$(EXE) \
+ $(bindir)/gpart$(EXE) $(bindir)/gscat$(EXE) \
+ $(bindir)/gtst$(EXE) $(bindir)/mcv$(EXE) \
+ $(bindir)/mmk_m2$(EXE) $(bindir)/mmk_m3$(EXE) \
+ $(bindir)/mord$(EXE) $(bindir)/mtst$(EXE)
+
+
ptinstall :
-$(CP) dggath$(EXE) dgmap$(EXE) dgord$(EXE) dgscat$(EXE) dgtst$(EXE) $(bindir)
-$(RM) $(bindir)/dgpart$(EXE)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -101,8 +101,8 @@
! elastic-poroelastic coupling surface
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: coupling_el_po_normal
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: coupling_el_po_jacobian2Dw
- integer, dimension(:,:,:), allocatable :: coupling_el_po_ijk
- integer, dimension(:), allocatable :: coupling_el_po_ispec
+ integer, dimension(:,:,:), allocatable :: coupling_el_po_ijk,coupling_po_el_ijk
+ integer, dimension(:), allocatable :: coupling_el_po_ispec,coupling_po_el_ispec
integer :: num_coupling_el_po_faces
! Moho mesh
@@ -438,8 +438,8 @@
coupling_ac_po_ijk,coupling_ac_po_ispec, &
num_coupling_ac_po_faces, &
coupling_el_po_normal,coupling_el_po_jacobian2Dw, &
- coupling_el_po_ijk,coupling_el_po_ispec, &
- num_coupling_el_po_faces, &
+ coupling_el_po_ijk,coupling_po_el_ijk,coupling_el_po_ispec, &
+ coupling_po_el_ispec,num_coupling_el_po_faces, &
num_interfaces_ext_mesh,my_neighbours_ext_mesh,nibool_interfaces_ext_mesh, &
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
prname,SAVE_MESH_FILES, &
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -425,7 +425,6 @@
integer, dimension(:), allocatable :: poroelastic_flag,acoustic_flag,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 = &
@@ -474,8 +473,6 @@
if( ier /= 0 ) stop 'error allocating array acoustic_flag'
allocate(test_flag(nglob),stat=ier)
if( ier /= 0 ) stop 'error allocating array test_flag'
- allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocating array mask_ibool'
poroelastic_flag(:) = 0
acoustic_flag(:) = 0
test_flag(:) = 0
@@ -538,10 +535,11 @@
! loops over all element faces and
! counts number of coupling faces between acoustic and poroelastic elements
- mask_ibool(:) = .false.
inum = 0
do ispec=1,nspec
+ if(ispec_is_poroelastic(ispec)) then
+
! loops over each face
do iface_ref= 1, 6
@@ -564,31 +562,10 @@
acoustic_flag( iglob_corners_ref(2) ) >= 1 .and. &
acoustic_flag( iglob_corners_ref(3) ) >= 1 .and. &
acoustic_flag( iglob_corners_ref(4) ) >= 1) then
- ! checks if face is has an poroelastic side
- if( poroelastic_flag( iglob_corners_ref(1) ) >= 1 .and. &
- poroelastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
- poroelastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
- poroelastic_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) .and. &
- ( acoustic_flag(iglob_midpoint) >= 1 .and. poroelastic_flag(iglob_midpoint) >= 1) ) then
-
- ! gets face GLL points i,j,k indices from element face
+ ! gets face GLL points i,j,k indices from poroelastic 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, &
@@ -605,9 +582,8 @@
ibool,nspec,nglob, &
xstore_dummy,ystore_dummy,zstore_dummy, &
normal_face(:,i,j) )
- ! makes sure that it always points away from acoustic element,
- ! otherwise switch direction
- if( ispec_is_poroelastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
+ ! reverse the sign, we know we are in a poroelastic element
+ normal_face(:,i,j) = - normal_face(:,i,j)
enddo
enddo
@@ -620,37 +596,24 @@
! 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?
+ ! we need to store local i,j,k,ispec info
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 ! poroelastic_flag
endif ! acoustic_flag
enddo ! iface_ref
+ endif ! ispec_is_poroelastic
enddo ! ispec
! stores completed coupling face informations
!
-! note: no need to store material parameters on these coupling points
-! for acoustic-poroelastic interface
+! note: for this coupling we need to have access to porous properties. The construction is such
+! that i,j,k, & face correspond to poroelastic interface. Note that the normal
+! is pointing outward the acoustic element
num_coupling_ac_po_faces = inum
allocate(coupling_ac_po_normal(NDIM,NGLLSQUARE,num_coupling_ac_po_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array coupling_ac_po_normal'
@@ -711,20 +674,19 @@
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 :: ijk_face_po(3,NGLLX,NGLLY), ijk_face_el(3,NGLLX,NGLLY)
+ integer,dimension(:,:,:),allocatable :: tmp_ijk,tmp_ijk_el
+ integer,dimension(:),allocatable :: tmp_ispec,tmp_ispec_el
- 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,dimension(NGNOD2D) :: iglob_corners_ref,iglob_corners_ref_el
+ integer :: ispec,i,j,k,igll,ier,iglob,ispec_el,ispec_ref_el
+ integer :: inum,iface_ref,iface_ref_el,iface_el,icorner,iglob_midpoint ! iface,ispec_neighbor
integer :: count_poroelastic,count_elastic
! mpi interface communication
integer, dimension(:), allocatable :: poroelastic_flag,elastic_flag,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 = &
@@ -759,10 +721,16 @@
if( ier /= 0 ) stop 'error allocating array tmp_jacobian2Dw'
allocate(tmp_ijk(3,NGLLSQUARE,nspec*6),stat=ier)
if( ier /= 0 ) stop 'error allocating array tmp_ijk'
+ allocate(tmp_ijk_el(3,NGLLSQUARE,nspec*6),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_ijk_el'
allocate(tmp_ispec(nspec*6),stat=ier)
if( ier /= 0 ) stop 'error allocating array tmp_ispec'
+ allocate(tmp_ispec_el(nspec*6),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_ispec_el'
tmp_ispec(:) = 0
+ tmp_ispec_el(:) = 0
tmp_ijk(:,:,:) = 0
+ tmp_ijk_el(:,:,:) = 0
tmp_normal(:,:,:) = 0.0
tmp_jacobian2Dw(:,:) = 0.0
@@ -773,8 +741,6 @@
if( ier /= 0 ) stop 'error allocating array elastic_flag'
allocate(test_flag(nglob),stat=ier)
if( ier /= 0 ) stop 'error allocating array test_flag'
- allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocating array mask_ibool'
poroelastic_flag(:) = 0
elastic_flag(:) = 0
test_flag(:) = 0
@@ -837,10 +803,11 @@
! loops over all element faces and
! counts number of coupling faces between elastic and poroelastic elements
- mask_ibool(:) = .false.
inum = 0
do ispec=1,nspec
+ if(ispec_is_poroelastic(ispec)) then
+
! loops over each face
do iface_ref= 1, 6
@@ -856,6 +823,7 @@
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 elastic side
@@ -863,30 +831,32 @@
elastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
elastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
elastic_flag( iglob_corners_ref(4) ) >= 1) then
- ! checks if face is has an poroelastic side
- if( poroelastic_flag( iglob_corners_ref(1) ) >= 1 .and. &
- poroelastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
- poroelastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
- poroelastic_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)
+ ! need to find elastic element for coupling
+ do ispec_el=1,nspec
+ if(ispec_is_elastic(ispec_el))then
+ do iface_el=6,1,-1
+ ! takes indices of corners of reference face
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface_el)
+ j = iface_all_corner_ijk(2,icorner,iface_el)
+ k = iface_all_corner_ijk(3,icorner,iface_el)
+ ! global reference indices
+ iglob_corners_ref_el(icorner) = ibool(i,j,k,ispec_el)
- ! checks if points on this face are masked already
- if( .not. mask_ibool(iglob_midpoint) .and. &
- ( elastic_flag(iglob_midpoint) >= 1 .and. poroelastic_flag(iglob_midpoint) >= 1) ) then
+ enddo
- ! gets face GLL points i,j,k indices from element face
- call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
+ if ( (iglob_corners_ref(1) == iglob_corners_ref_el(3)) .and. &
+ (iglob_corners_ref(3) == iglob_corners_ref_el(1)) ) then
- ! 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
+ iface_ref_el = iface_el ![CM]: for some reason this shows a wrong orientation
+ ! but the calcul is ok.
+ ispec_ref_el = ispec_el
+
+ ! gets face GLL points i,j,k indices from poroelastic element face
+ call get_element_face_gll_indices(iface_ref,ijk_face_po,NGLLX,NGLLY)
+ ! gets face GLL points i,j,k indices from elastic element face
+ call get_element_face_gll_indices(iface_ref_el,ijk_face_el,NGLLX,NGLLY)
! gets face GLL 2Djacobian, weighted from element face
call get_jacobian_boundary_face(myrank,nspec, &
@@ -895,61 +865,51 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
- ! normal convention: points away from elastic, reference element
- ! switch normal direction if necessary
+ ! normal convention: points away from poroelastic, reference element
do j=1,NGLLY
do i=1,NGLLX
- ! directs normals such that they point outwards of element
+ ! directs normals such that they point outwards of poroelastic 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 elastic element,
- ! otherwise switch direction
- if( ispec_is_poroelastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
enddo
enddo
! stores informations about this face
inum = inum + 1
tmp_ispec(inum) = ispec
+ tmp_ispec_el(inum) = ispec_ref_el
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)
+ ! we need to store local i,j,k,ispec info
+ tmp_ijk(:,igll,inum) = ijk_face_po(:,i,j)
+ tmp_ijk_el(:,igll,inum) = ijk_face_el(:,NGLLY-j+1,NGLLX-i+1)
! 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 ! poroelastic_flag
+ endif ! if
+
+ enddo ! do iface_ref_el=1,6
+ endif ! if(ispec_is_elastic(ispec_el))then
+ enddo ! do ispec_el=1,nspec
endif ! elastic_flag
enddo ! iface_ref
+ endif ! ispec_is_poroelastic
enddo ! ispec
! stores completed coupling face informations
!
-! note: no need to store material parameters on these coupling points
-! for elastic-poroelastic interface
+! note: for this coupling we need to have access to porous properties. The construction is such
+! that i,j,k, & face correspond to poroelastic interface. Note that the normal
+! is pointing outward the poroelastic element
num_coupling_el_po_faces = inum
allocate(coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array coupling_el_po_normal'
@@ -957,13 +917,19 @@
if( ier /= 0 ) stop 'error allocating array coupling_el_po_jacobian2Dw'
allocate(coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array coupling_el_po_ijk'
+ allocate(coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_po_el_ijk'
allocate(coupling_el_po_ispec(num_coupling_el_po_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array coupling_el_po_ispec'
+ allocate(coupling_po_el_ispec(num_coupling_el_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_po_el_ispec'
do inum = 1,num_coupling_el_po_faces
coupling_el_po_normal(:,:,inum) = tmp_normal(:,:,inum)
coupling_el_po_jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
coupling_el_po_ijk(:,:,inum) = tmp_ijk(:,:,inum)
+ coupling_po_el_ijk(:,:,inum) = tmp_ijk_el(:,:,inum)
coupling_el_po_ispec(inum) = tmp_ispec(inum)
+ coupling_po_el_ispec(inum) = tmp_ispec_el(inum)
enddo
! user output
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -51,8 +51,8 @@
coupling_ac_po_ijk,coupling_ac_po_ispec, &
num_coupling_ac_po_faces, &
coupling_el_po_normal,coupling_el_po_jacobian2Dw, &
- coupling_el_po_ijk,coupling_el_po_ispec, &
- num_coupling_el_po_faces, &
+ coupling_el_po_ijk,coupling_po_el_ijk,coupling_el_po_ispec, &
+ coupling_po_el_ispec,num_coupling_el_po_faces, &
num_interfaces_ext_mesh,my_neighbours_ext_mesh,nibool_interfaces_ext_mesh, &
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
prname,SAVE_MESH_FILES, &
@@ -133,7 +133,9 @@
real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
real(kind=CUSTOM_REAL) :: coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces)
integer :: coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
+ integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
! MPI interfaces
integer :: num_interfaces_ext_mesh
@@ -298,7 +300,9 @@
write(IOUT) num_coupling_el_po_faces
if( num_coupling_el_po_faces > 0 ) then
write(IOUT) coupling_el_po_ispec
+ write(IOUT) coupling_po_el_ispec
write(IOUT) coupling_el_po_ijk
+ write(IOUT) coupling_po_el_ijk
write(IOUT) coupling_el_po_jacobian2Dw
write(IOUT) coupling_el_po_normal
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -313,7 +313,7 @@
!
subroutine check_mesh_resolution_poro(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
- kappastore,mustore,rho_vp,rho_vs, &
+ kappastore,mustore, &
DT, model_speed_max,min_resolved_period, &
phistore,tortstore,rhoarraystore, &
rho_vpI,rho_vpII,rho_vsI )
@@ -327,7 +327,7 @@
include "constants.h"
integer :: NSPEC_AB,NGLOB_AB
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappastore,mustore,rho_vp,rho_vs
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappastore,mustore
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vpI,rho_vpII,rho_vsI
real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: phistore,tortstore
@@ -400,7 +400,7 @@
! determines minimum/maximum velocities within this element
call get_vpvs_minmax_poro(vpmin,vpmax,vp2min,vp2max,vsmin,vsmax,ispec,has_vs_zero, &
- has_vp2_zero,NSPEC_AB,kappastore,mustore,rho_vp,rho_vs, &
+ has_vp2_zero,NSPEC_AB,kappastore,mustore, &
phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI)
! min/max for whole cpu partition
@@ -459,7 +459,7 @@
if (has_vp2_zero) then
pmax = avg_distance / min( vpmin,vsmin ) * NPTS_PER_WAVELENGTH
else
- pmax = avg_distance / min( vpmin,vp2max,vsmin ) * NPTS_PER_WAVELENGTH
+ pmax = avg_distance / min( vpmin,vp2min,vsmin ) * NPTS_PER_WAVELENGTH
endif
pmax_glob = max(pmax_glob,pmax)
@@ -697,11 +697,12 @@
subroutine get_vpvs_minmax_poro(vpmin,vpmax,vp2min,vp2max,vsmin,vsmax,ispec,has_vs_zero, &
- has_vp2_zero,NSPEC_AB,kappastore,mustore,rho_vp,rho_vs, &
+ has_vp2_zero,NSPEC_AB,kappastore,mustore, &
phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI)
! calculates the min/max size of the specified element (ispec) for poroelastic domains
-
+! [CM] Note: in case of coupled acoustic-poro-elastic, rho_vpI,rho_vpII,rho_vsI, etc have
+! been appropriately defined in /src/generate_databases/get_model.f90
implicit none
include "constants.h"
@@ -714,7 +715,7 @@
integer :: NSPEC_AB
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
- kappastore,mustore,rho_vp,rho_vs
+ kappastore,mustore
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vpI,rho_vpII,rho_vsI
real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: phistore,tortstore
@@ -732,10 +733,7 @@
vsmax = -HUGEVAL
! vp
- where( rho_vp(:,:,:,ispec) > TINYVAL )
- vp_elem(:,:,:) = (FOUR_THIRDS * mustore(:,:,:,ispec) &
- + kappastore(:,:,:,ispec)) / rho_vp(:,:,:,ispec)
- elsewhere(rho_vpI(:,:,:,ispec) > TINYVAL)
+ where(rho_vpI(:,:,:,ispec) > TINYVAL)
vp_elem(:,:,:) = rho_vpI(:,:,:,ispec)/( ((1._CUSTOM_REAL - phistore(:,:,:,ispec))* &
rhoarraystore(1,:,:,:,ispec) + phistore(:,:,:,ispec)*rhoarraystore(2,:,:,:,ispec)) - &
phistore(:,:,:,ispec)/tortstore(:,:,:,ispec) * rhoarraystore(2,:,:,:,ispec) )
@@ -770,9 +768,7 @@
vp2max = max(vp2max,val_max(1))
! vs
- where( rho_vs(:,:,:,ispec) > TINYVAL )
- vs_elem(:,:,:) = mustore(:,:,:,ispec) / rho_vs(:,:,:,ispec)
- elsewhere(rho_vsI(:,:,:,ispec) > TINYVAL)
+ where(rho_vsI(:,:,:,ispec) > TINYVAL)
vs_elem(:,:,:) = rho_vsI(:,:,:,ispec) / ( ((1._CUSTOM_REAL - phistore(:,:,:,ispec))* &
rhoarraystore(1,:,:,:,ispec) + phistore(:,:,:,ispec)*rhoarraystore(2,:,:,:,ispec)) - &
phistore(:,:,:,ispec)/tortstore(:,:,:,ispec) * rhoarraystore(2,:,:,:,ispec) )
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2011-11-29 04:02:32 UTC (rev 19245)
@@ -123,8 +123,10 @@
$O/compute_add_sources_poroelastic.o \
$O/compute_coupling_acoustic_el.o \
$O/compute_coupling_elastic_ac.o \
+ $O/compute_coupling_elastic_po.o \
$O/compute_coupling_acoustic_po.o \
$O/compute_coupling_poroelastic_ac.o \
+ $O/compute_coupling_poroelastic_el.o \
$O/compute_stacey_acoustic.o \
$O/compute_stacey_elastic.o \
$O/compute_gradient.o \
@@ -314,12 +316,18 @@
$O/compute_coupling_elastic_ac.o: $(SHARED)constants.h compute_coupling_elastic_ac.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_elastic_ac.o compute_coupling_elastic_ac.f90
+$O/compute_coupling_elastic_po.o: $(SHARED)constants.h compute_coupling_elastic_po.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_elastic_po.o compute_coupling_elastic_po.f90
+
$O/compute_coupling_acoustic_po.o: $(SHARED)constants.h compute_coupling_acoustic_po.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_acoustic_po.o compute_coupling_acoustic_po.f90
$O/compute_coupling_poroelastic_ac.o: $(SHARED)constants.h compute_coupling_poroelastic_ac.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_poroelastic_ac.o compute_coupling_poroelastic_ac.f90
+$O/compute_coupling_poroelastic_el.o: $(SHARED)constants.h compute_coupling_poroelastic_el.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling_poroelastic_el.o compute_coupling_poroelastic_el.f90
+
$O/compute_stacey_acoustic.o: $(SHARED)constants.h compute_stacey_acoustic.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_stacey_acoustic.o compute_stacey_acoustic.f90
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -0,0 +1,496 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+! for elastic solver
+
+ subroutine compute_coupling_elastic_po(NSPEC_AB,NGLOB_AB,ibool,&
+ displs_poroelastic,accels_poroelastic,displw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,&
+ displ,accel,kappastore, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store, &
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_coupling_el_po_faces, &
+ coupling_el_po_ispec,coupling_po_el_ispec, &
+ coupling_el_po_ijk,coupling_po_el_ijk, &
+ coupling_el_po_normal, &
+ coupling_el_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+
+! returns the updated accelerations array: accel
+
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacements, etc
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,accels_poroelastic,&
+ displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+
+! global indexing
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ integer :: SIMULATION_TYPE
+ integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+
+! elastic-poroelastic coupling surface
+ integer :: num_coupling_el_po_faces
+ real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
+ integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
+
+! properties
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ etastore,phistore,tortstore,jacobian
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappaarraystore
+ real(kind=CUSTOM_REAL), dimension(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: permstore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ kappastore,mustore
+
+! anisotropy
+ logical :: ANISOTROPY
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! local parameters
+ real(kind=CUSTOM_REAL) :: sigmap,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,phil,tortl,rhol_bar
+ real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) :: kappal_s
+ real(kind=CUSTOM_REAL) :: etal_f,kappal_f
+ real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot
+ real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
+
+! local anisotropy parameters
+ real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+ integer :: iface,igll,ispec_po,ispec_el,iglob,iglob_el,iglob_po
+ integer :: i,j,k,l
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+ real(kind=CUSTOM_REAL) tempx1ls,tempx2ls,tempx3ls,tempx1lw,tempx2lw,tempx3lw
+ real(kind=CUSTOM_REAL) tempy1ls,tempy2ls,tempy3ls,tempy1lw,tempy2lw,tempy3lw
+ real(kind=CUSTOM_REAL) tempz1ls,tempz2ls,tempz3ls,tempz1lw,tempz2lw,tempz3lw
+
+ real(kind=CUSTOM_REAL) :: duxdxl,duydxl,duzdxl,duxdyl,duydyl,duzdyl,duxdzl,duydzl,duzdzl
+ real(kind=CUSTOM_REAL) :: dwxdxl,dwydxl,dwzdxl,dwxdyl,dwydyl,dwzdyl,dwxdzl,dwydzl,dwzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl_plus_duzdzl,dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+
+! Jacobian matrix and determinant
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+
+! loops on all coupling faces
+ do iface = 1,num_coupling_el_po_faces
+
+ ! gets corresponding poro/elastic spectral element
+ ispec_po = coupling_el_po_ispec(iface)
+ ispec_el = coupling_po_el_ispec(iface)
+
+ if( ispec_is_inner(ispec_el) .eqv. phase_is_inner ) then
+
+ ! loops over common GLL points
+ do igll = 1, NGLLSQUARE
+
+ !-----------------------
+ ! from the poroelastic side
+ !-----------------------
+ i = coupling_el_po_ijk(1,igll,iface)
+ j = coupling_el_po_ijk(2,igll,iface)
+ k = coupling_el_po_ijk(3,igll,iface)
+
+ iglob_po = ibool(i,j,k,ispec_po)
+
+! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec_po)
+ tortl = tortstore(i,j,k,ispec_po)
+!solid properties
+ kappal_s = kappaarraystore(1,i,j,k,ispec_po)
+ rhol_s = rhoarraystore(1,i,j,k,ispec_po)
+!fluid properties
+ kappal_f = kappaarraystore(2,i,j,k,ispec_po)
+ rhol_f = rhoarraystore(2,i,j,k,ispec_po)
+!frame properties
+ mul_fr = mustore(i,j,k,ispec_po)
+ kappal_fr = kappaarraystore(3,i,j,k,ispec_po)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+!Biot coefficients for the input phi
+ D_biot = kappal_s*(1._CUSTOM_REAL + phil*(kappal_s/kappal_f - 1._CUSTOM_REAL))
+ H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + &
+ kappal_fr + 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
+ C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
+ M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
+
+!T = G:grad u_s + C_biot div w I
+!and T_f = C_biot div u_s I + M_biot div w I
+ mul_G = mul_fr
+ lambdal_G = H_biot - 2._CUSTOM_REAL*mul_fr
+ lambdalplus2mul_G = lambdal_G + 2._CUSTOM_REAL*mul_G
+
+! derivative along x,y,z for u_s and w
+ tempx1ls = 0.
+ tempx2ls = 0.
+ tempx3ls = 0.
+
+ tempy1ls = 0.
+ tempy2ls = 0.
+ tempy3ls = 0.
+
+ tempz1ls = 0.
+ tempz2ls = 0.
+ tempz3ls = 0.
+
+ tempx1lw = 0.
+ tempx2lw = 0.
+ tempx3lw = 0.
+
+ tempy1lw = 0.
+ tempy2lw = 0.
+ tempy3lw = 0.
+
+ tempz1lw = 0.
+ tempz2lw = 0.
+ tempz3lw = 0.
+
+! first double loop over GLL points to compute and store gradients
+ do l = 1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec_po)
+ tempx1ls = tempx1ls + displs_poroelastic(1,iglob)*hp1
+ tempy1ls = tempy1ls + displs_poroelastic(2,iglob)*hp1
+ tempz1ls = tempz1ls + displs_poroelastic(3,iglob)*hp1
+ tempx1lw = tempx1lw + displw_poroelastic(1,iglob)*hp1
+ tempy1lw = tempy1lw + displw_poroelastic(2,iglob)*hp1
+ tempz1lw = tempz1lw + displw_poroelastic(3,iglob)*hp1
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do
+ !l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec_po)
+ tempx2ls = tempx2ls + displs_poroelastic(1,iglob)*hp2
+ tempy2ls = tempy2ls + displs_poroelastic(2,iglob)*hp2
+ tempz2ls = tempz2ls + displs_poroelastic(3,iglob)*hp2
+ tempx2lw = tempx2lw + displw_poroelastic(1,iglob)*hp2
+ tempy2lw = tempy2lw + displw_poroelastic(2,iglob)*hp2
+ tempz2lw = tempz2lw + displw_poroelastic(3,iglob)*hp2
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do
+ !l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec_po)
+ tempx3ls = tempx3ls + displs_poroelastic(1,iglob)*hp3
+ tempy3ls = tempy3ls + displs_poroelastic(2,iglob)*hp3
+ tempz3ls = tempz3ls + displs_poroelastic(3,iglob)*hp3
+ tempx3lw = tempx3lw + displw_poroelastic(1,iglob)*hp3
+ tempy3lw = tempy3lw + displw_poroelastic(2,iglob)*hp3
+ tempz3lw = tempz3lw + displw_poroelastic(3,iglob)*hp3
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ enddo
+
+ xixl = xix(i,j,k,ispec_po)
+ xiyl = xiy(i,j,k,ispec_po)
+ xizl = xiz(i,j,k,ispec_po)
+ etaxl = etax(i,j,k,ispec_po)
+ etayl = etay(i,j,k,ispec_po)
+ etazl = etaz(i,j,k,ispec_po)
+ gammaxl = gammax(i,j,k,ispec_po)
+ gammayl = gammay(i,j,k,ispec_po)
+ gammazl = gammaz(i,j,k,ispec_po)
+ jacobianl = jacobian(i,j,k,ispec_po)
+
+! derivatives of displacement
+ duxdxl = xixl*tempx1ls + etaxl*tempx2ls + gammaxl*tempx3ls
+ duxdyl = xiyl*tempx1ls + etayl*tempx2ls + gammayl*tempx3ls
+ duxdzl = xizl*tempx1ls + etazl*tempx2ls + gammazl*tempx3ls
+
+ duydxl = xixl*tempy1ls + etaxl*tempy2ls + gammaxl*tempy3ls
+ duydyl = xiyl*tempy1ls + etayl*tempy2ls + gammayl*tempy3ls
+ duydzl = xizl*tempy1ls + etazl*tempy2ls + gammazl*tempy3ls
+
+ duzdxl = xixl*tempz1ls + etaxl*tempz2ls + gammaxl*tempz3ls
+ duzdyl = xiyl*tempz1ls + etayl*tempz2ls + gammayl*tempz3ls
+ duzdzl = xizl*tempz1ls + etazl*tempz2ls + gammazl*tempz3ls
+
+ dwxdxl = xixl*tempx1lw + etaxl*tempx2lw + gammaxl*tempx3lw
+ dwxdyl = xiyl*tempx1lw + etayl*tempx2lw + gammayl*tempx3lw
+ dwxdzl = xizl*tempx1lw + etazl*tempx2lw + gammazl*tempx3lw
+
+ dwydxl = xixl*tempy1lw + etaxl*tempy2lw + gammaxl*tempy3lw
+ dwydyl = xiyl*tempy1lw + etayl*tempy2lw + gammayl*tempy3lw
+ dwydzl = xizl*tempy1lw + etazl*tempy2lw + gammazl*tempy3lw
+
+ dwzdxl = xixl*tempz1lw + etaxl*tempz2lw + gammaxl*tempz3lw
+ dwzdyl = xiyl*tempz1lw + etayl*tempz2lw + gammayl*tempz3lw
+ dwzdzl = xizl*tempz1lw + etazl*tempz2lw + gammazl*tempz3lw
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl_plus_duzdzl = duxdxl + duydyl + duzdzl
+ dwxdxl_plus_dwydyl_plus_dwzdzl = dwxdxl + dwydyl + dwzdzl
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute stress tensor (include attenuation or anisotropy if needed)
+
+! if(VISCOATTENUATION) then
+!chris:check
+
+! Dissipation only controlled by frame share attenuation in poroelastic (see
+! Morency & Tromp, GJI 2008).
+! attenuation is implemented following the memory variable formulation of
+! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
+! vol. 58(1), p. 110-120 (1993). More details can be found in
+! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a
+! linear
+! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611
+! (1988).
+
+! else
+
+! no attenuation
+ sigma_xx = lambdalplus2mul_G*duxdxl + lambdal_G*duydyl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_yy = lambdalplus2mul_G*duydyl + lambdal_G*duxdxl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_zz = lambdalplus2mul_G*duzdzl + lambdal_G*duxdxl_plus_duydyl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ sigma_xy = mul_G*duxdyl_plus_duydxl
+ sigma_xz = mul_G*duzdxl_plus_duxdzl
+ sigma_yz = mul_G*duzdyl_plus_duydzl
+
+
+ !-----------------------
+ ! from the elastic side
+ !-----------------------
+ i = coupling_po_el_ijk(1,igll,iface)
+ j = coupling_po_el_ijk(2,igll,iface)
+ k = coupling_po_el_ijk(3,igll,iface)
+
+ ! gets global index of this common GLL point
+ ! (note: should be the same as for corresponding
+ ! i',j',k',ispec_poroelastic or ispec_elastic )
+ iglob_el = ibool(i,j,k,ispec_el)
+ if (iglob_el .ne. iglob_po) stop 'poroelastic-elastic coupling error'
+ tempx1l = 0.
+ tempx2l = 0.
+ tempx3l = 0.
+
+ tempy1l = 0.
+ tempy2l = 0.
+ tempy3l = 0.
+
+ tempz1l = 0.
+ tempz2l = 0.
+ tempz3l = 0.
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec_el)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec_el)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec_el)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
+
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec_el)
+ xiyl = xiy(i,j,k,ispec_el)
+ xizl = xiz(i,j,k,ispec_el)
+ etaxl = etax(i,j,k,ispec_el)
+ etayl = etay(i,j,k,ispec_el)
+ etazl = etaz(i,j,k,ispec_el)
+ gammaxl = gammax(i,j,k,ispec_el)
+ gammayl = gammay(i,j,k,ispec_el)
+ gammazl = gammaz(i,j,k,ispec_el)
+ jacobianl = jacobian(i,j,k,ispec_el)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ kappal = kappastore(i,j,k,ispec_el)
+ mul = mustore(i,j,k,ispec_el)
+
+ ! full anisotropic case, stress calculations
+ if(ANISOTROPY) then
+ c11 = c11store(i,j,k,ispec_el)
+ c12 = c12store(i,j,k,ispec_el)
+ c13 = c13store(i,j,k,ispec_el)
+ c14 = c14store(i,j,k,ispec_el)
+ c15 = c15store(i,j,k,ispec_el)
+ c16 = c16store(i,j,k,ispec_el)
+ c22 = c22store(i,j,k,ispec_el)
+ c23 = c23store(i,j,k,ispec_el)
+ c24 = c24store(i,j,k,ispec_el)
+ c25 = c25store(i,j,k,ispec_el)
+ c26 = c26store(i,j,k,ispec_el)
+ c33 = c33store(i,j,k,ispec_el)
+ c34 = c34store(i,j,k,ispec_el)
+ c35 = c35store(i,j,k,ispec_el)
+ c36 = c36store(i,j,k,ispec_el)
+ c44 = c44store(i,j,k,ispec_el)
+ c45 = c45store(i,j,k,ispec_el)
+ c46 = c46store(i,j,k,ispec_el)
+ c55 = c55store(i,j,k,ispec_el)
+ c56 = c56store(i,j,k,ispec_el)
+ c66 = c66store(i,j,k,ispec_el)
+
+ sigma_xx = sigma_xx + c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+ sigma_yy = sigma_yy + c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+ sigma_zz = sigma_zz + c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+ sigma_xy = sigma_xy + c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+ sigma_xz = sigma_xz + c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+ sigma_yz = sigma_yz + c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+
+ else
+
+ ! isotropic case
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+ ! compute stress sigma
+ sigma_xx = sigma_xx + lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = sigma_yy + lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = sigma_zz + lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = sigma_xy + mul*duxdyl_plus_duydxl
+ sigma_xz = sigma_xz + mul*duzdxl_plus_duxdzl
+ sigma_yz = sigma_yz + mul*duzdyl_plus_duydzl
+
+ endif ! ANISOTROPY
+
+ ! gets associated normal on GLL point
+ ! (note convention: pointing outwards of poroelastic element)
+ nx = coupling_el_po_normal(1,igll,iface)
+ ny = coupling_el_po_normal(2,igll,iface)
+ nz = coupling_el_po_normal(3,igll,iface)
+
+ ! gets associated, weighted 2D jacobian
+ ! (note: should be the same for poroelastic and elastic element)
+ jacobianw = coupling_el_po_jacobian2Dw(igll,iface)
+
+ ! continuity of displacement and traction on global point
+ !
+ ! note: continuity of displacement is enforced after the velocity update
+ accel(1,iglob_el) = accel(1,iglob_el) - jacobianw* &
+ ( sigma_xx*nx + sigma_xy*ny + sigma_xz*nz )/2.d0
+
+ accel(2,iglob_el) = accel(2,iglob_el) - jacobianw* &
+ ( sigma_xy*nx + sigma_yy*ny + sigma_yz*nz )/2.d0
+
+ accel(3,iglob_el) = accel(3,iglob_el) - jacobianw* &
+ ( sigma_xz*nx + sigma_yz*ny + sigma_zz*nz )/2.d0
+
+ enddo ! igll
+
+ endif
+
+ enddo ! iface
+
+end subroutine compute_coupling_elastic_po
+
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -0,0 +1,505 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+! for poroelastic solver
+
+ subroutine compute_coupling_poroelastic_el(NSPEC_AB,NGLOB_AB,ibool,&
+ displs_poroelastic,accels_poroelastic,displw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,&
+ displ,accel,kappastore, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store, &
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_coupling_el_po_faces, &
+ coupling_el_po_ispec,coupling_po_el_ispec, &
+ coupling_el_po_ijk,coupling_po_el_ijk, &
+ coupling_el_po_normal, &
+ coupling_el_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+
+! returns the updated accelerations array: accels_poroelatsic (& accelw_poroelastic )
+
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacements, etc
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,accels_poroelastic,&
+ displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+
+! global indexing
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ integer :: SIMULATION_TYPE
+ integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+
+! elastic-poroelastic coupling surface
+ integer :: num_coupling_el_po_faces
+ real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
+ integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
+
+! properties
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ etastore,phistore,tortstore,jacobian
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappaarraystore
+ real(kind=CUSTOM_REAL), dimension(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: permstore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ kappastore,mustore
+
+! anisotropy
+ logical :: ANISOTROPY
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! local parameters
+ real(kind=CUSTOM_REAL) :: sigmap,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,phil,tortl,rhol_bar
+ real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL) :: kappal_s
+ real(kind=CUSTOM_REAL) :: etal_f,kappal_f
+ real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot
+ real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
+
+! local anisotropy parameters
+ real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+ integer :: iface,igll,ispec_po,ispec_el,iglob,iglob_po,iglob_el
+ integer :: i,j,k,l
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+ real(kind=CUSTOM_REAL) tempx1ls,tempx2ls,tempx3ls,tempx1lw,tempx2lw,tempx3lw
+ real(kind=CUSTOM_REAL) tempy1ls,tempy2ls,tempy3ls,tempy1lw,tempy2lw,tempy3lw
+ real(kind=CUSTOM_REAL) tempz1ls,tempz2ls,tempz3ls,tempz1lw,tempz2lw,tempz3lw
+
+ real(kind=CUSTOM_REAL) :: duxdxl,duydxl,duzdxl,duxdyl,duydyl,duzdyl,duxdzl,duydzl,duzdzl
+ real(kind=CUSTOM_REAL) :: dwxdxl,dwydxl,dwzdxl,dwxdyl,dwydyl,dwzdyl,dwxdzl,dwydzl,dwzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl_plus_duzdzl,dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+
+! Jacobian matrix and determinant
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+
+! loops on all coupling faces
+ do iface = 1,num_coupling_el_po_faces
+
+ ! gets corresponding poro/elastic spectral element
+ ispec_po = coupling_el_po_ispec(iface)
+ ispec_el = coupling_po_el_ispec(iface)
+
+ if( ispec_is_inner(ispec_po) .eqv. phase_is_inner ) then
+
+ ! loops over common GLL points
+ do igll = 1, NGLLSQUARE
+
+ !-----------------------
+ ! from the elastic side
+ !-----------------------
+ i = coupling_po_el_ijk(1,igll,iface)
+ j = coupling_po_el_ijk(2,igll,iface)
+ k = coupling_po_el_ijk(3,igll,iface)
+
+ ! gets global index of this common GLL point
+ ! (note: should be the same as for corresponding i',j',k',ispec_poroelastic or ispec_elastic )
+ iglob_el = ibool(i,j,k,ispec_el)
+
+ tempx1l = 0.
+ tempx2l = 0.
+ tempx3l = 0.
+
+ tempy1l = 0.
+ tempy2l = 0.
+ tempy3l = 0.
+
+ tempz1l = 0.
+ tempz2l = 0.
+ tempz3l = 0.
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec_el)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec_el)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec_el)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
+
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec_el)
+ xiyl = xiy(i,j,k,ispec_el)
+ xizl = xiz(i,j,k,ispec_el)
+ etaxl = etax(i,j,k,ispec_el)
+ etayl = etay(i,j,k,ispec_el)
+ etazl = etaz(i,j,k,ispec_el)
+ gammaxl = gammax(i,j,k,ispec_el)
+ gammayl = gammay(i,j,k,ispec_el)
+ gammazl = gammaz(i,j,k,ispec_el)
+ jacobianl = jacobian(i,j,k,ispec_el)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ kappal = kappastore(i,j,k,ispec_el)
+ mul = mustore(i,j,k,ispec_el)
+
+ ! full anisotropic case, stress calculations
+ if(ANISOTROPY) then
+ c11 = c11store(i,j,k,ispec_el)
+ c12 = c12store(i,j,k,ispec_el)
+ c13 = c13store(i,j,k,ispec_el)
+ c14 = c14store(i,j,k,ispec_el)
+ c15 = c15store(i,j,k,ispec_el)
+ c16 = c16store(i,j,k,ispec_el)
+ c22 = c22store(i,j,k,ispec_el)
+ c23 = c23store(i,j,k,ispec_el)
+ c24 = c24store(i,j,k,ispec_el)
+ c25 = c25store(i,j,k,ispec_el)
+ c26 = c26store(i,j,k,ispec_el)
+ c33 = c33store(i,j,k,ispec_el)
+ c34 = c34store(i,j,k,ispec_el)
+ c35 = c35store(i,j,k,ispec_el)
+ c36 = c36store(i,j,k,ispec_el)
+ c44 = c44store(i,j,k,ispec_el)
+ c45 = c45store(i,j,k,ispec_el)
+ c46 = c46store(i,j,k,ispec_el)
+ c55 = c55store(i,j,k,ispec_el)
+ c56 = c56store(i,j,k,ispec_el)
+ c66 = c66store(i,j,k,ispec_el)
+
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+
+ else
+
+ ! isotropic case
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+ endif ! ANISOTROPY
+
+
+ !-----------------------
+ ! from the poroelastic side
+ !-----------------------
+ i = coupling_el_po_ijk(1,igll,iface)
+ j = coupling_el_po_ijk(2,igll,iface)
+ k = coupling_el_po_ijk(3,igll,iface)
+
+ iglob_po = ibool(i,j,k,ispec_po)
+ if (iglob_el .ne. iglob_po) stop 'poroelastic-elastic coupling error'
+
+! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec_po)
+ tortl = tortstore(i,j,k,ispec_po)
+!solid properties
+ kappal_s = kappaarraystore(1,i,j,k,ispec_po)
+ rhol_s = rhoarraystore(1,i,j,k,ispec_po)
+!fluid properties
+ kappal_f = kappaarraystore(2,i,j,k,ispec_po)
+ rhol_f = rhoarraystore(2,i,j,k,ispec_po)
+!frame properties
+ mul_fr = mustore(i,j,k,ispec_po)
+ kappal_fr = kappaarraystore(3,i,j,k,ispec_po)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+!Biot coefficients for the input phi
+ D_biot = kappal_s*(1._CUSTOM_REAL + phil*(kappal_s/kappal_f - 1._CUSTOM_REAL))
+ H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + &
+ kappal_fr + 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
+ C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
+ M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
+
+!T = G:grad u_s + C_biot div w I
+!and T_f = C_biot div u_s I + M_biot div w I
+ mul_G = mul_fr
+ lambdal_G = H_biot - 2._CUSTOM_REAL*mul_fr
+ lambdalplus2mul_G = lambdal_G + 2._CUSTOM_REAL*mul_G
+
+! derivative along x,y,z for u_s and w
+ tempx1ls = 0.
+ tempx2ls = 0.
+ tempx3ls = 0.
+
+ tempy1ls = 0.
+ tempy2ls = 0.
+ tempy3ls = 0.
+
+ tempz1ls = 0.
+ tempz2ls = 0.
+ tempz3ls = 0.
+
+ tempx1lw = 0.
+ tempx2lw = 0.
+ tempx3lw = 0.
+
+ tempy1lw = 0.
+ tempy2lw = 0.
+ tempy3lw = 0.
+
+ tempz1lw = 0.
+ tempz2lw = 0.
+ tempz3lw = 0.
+
+! first double loop over GLL points to compute and store gradients
+ do l = 1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec_po)
+ tempx1ls = tempx1ls + displs_poroelastic(1,iglob)*hp1
+ tempy1ls = tempy1ls + displs_poroelastic(2,iglob)*hp1
+ tempz1ls = tempz1ls + displs_poroelastic(3,iglob)*hp1
+ tempx1lw = tempx1lw + displw_poroelastic(1,iglob)*hp1
+ tempy1lw = tempy1lw + displw_poroelastic(2,iglob)*hp1
+ tempz1lw = tempz1lw + displw_poroelastic(3,iglob)*hp1
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do
+ !l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec_po)
+ tempx2ls = tempx2ls + displs_poroelastic(1,iglob)*hp2
+ tempy2ls = tempy2ls + displs_poroelastic(2,iglob)*hp2
+ tempz2ls = tempz2ls + displs_poroelastic(3,iglob)*hp2
+ tempx2lw = tempx2lw + displw_poroelastic(1,iglob)*hp2
+ tempy2lw = tempy2lw + displw_poroelastic(2,iglob)*hp2
+ tempz2lw = tempz2lw + displw_poroelastic(3,iglob)*hp2
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do
+ !l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec_po)
+ tempx3ls = tempx3ls + displs_poroelastic(1,iglob)*hp3
+ tempy3ls = tempy3ls + displs_poroelastic(2,iglob)*hp3
+ tempz3ls = tempz3ls + displs_poroelastic(3,iglob)*hp3
+ tempx3lw = tempx3lw + displw_poroelastic(1,iglob)*hp3
+ tempy3lw = tempy3lw + displw_poroelastic(2,iglob)*hp3
+ tempz3lw = tempz3lw + displw_poroelastic(3,iglob)*hp3
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ enddo
+
+ xixl = xix(i,j,k,ispec_po)
+ xiyl = xiy(i,j,k,ispec_po)
+ xizl = xiz(i,j,k,ispec_po)
+ etaxl = etax(i,j,k,ispec_po)
+ etayl = etay(i,j,k,ispec_po)
+ etazl = etaz(i,j,k,ispec_po)
+ gammaxl = gammax(i,j,k,ispec_po)
+ gammayl = gammay(i,j,k,ispec_po)
+ gammazl = gammaz(i,j,k,ispec_po)
+ jacobianl = jacobian(i,j,k,ispec_po)
+
+! derivatives of displacement
+ duxdxl = xixl*tempx1ls + etaxl*tempx2ls + gammaxl*tempx3ls
+ duxdyl = xiyl*tempx1ls + etayl*tempx2ls + gammayl*tempx3ls
+ duxdzl = xizl*tempx1ls + etazl*tempx2ls + gammazl*tempx3ls
+
+ duydxl = xixl*tempy1ls + etaxl*tempy2ls + gammaxl*tempy3ls
+ duydyl = xiyl*tempy1ls + etayl*tempy2ls + gammayl*tempy3ls
+ duydzl = xizl*tempy1ls + etazl*tempy2ls + gammazl*tempy3ls
+
+ duzdxl = xixl*tempz1ls + etaxl*tempz2ls + gammaxl*tempz3ls
+ duzdyl = xiyl*tempz1ls + etayl*tempz2ls + gammayl*tempz3ls
+ duzdzl = xizl*tempz1ls + etazl*tempz2ls + gammazl*tempz3ls
+
+ dwxdxl = xixl*tempx1lw + etaxl*tempx2lw + gammaxl*tempx3lw
+ dwxdyl = xiyl*tempx1lw + etayl*tempx2lw + gammayl*tempx3lw
+ dwxdzl = xizl*tempx1lw + etazl*tempx2lw + gammazl*tempx3lw
+
+ dwydxl = xixl*tempy1lw + etaxl*tempy2lw + gammaxl*tempy3lw
+ dwydyl = xiyl*tempy1lw + etayl*tempy2lw + gammayl*tempy3lw
+ dwydzl = xizl*tempy1lw + etazl*tempy2lw + gammazl*tempy3lw
+
+ dwzdxl = xixl*tempz1lw + etaxl*tempz2lw + gammaxl*tempz3lw
+ dwzdyl = xiyl*tempz1lw + etayl*tempz2lw + gammayl*tempz3lw
+ dwzdzl = xizl*tempz1lw + etazl*tempz2lw + gammazl*tempz3lw
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl_plus_duzdzl = duxdxl + duydyl + duzdzl
+ dwxdxl_plus_dwydyl_plus_dwzdzl = dwxdxl + dwydyl + dwzdzl
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute stress tensor (include attenuation or anisotropy if needed)
+
+! if(VISCOATTENUATION) then
+!chris:check
+
+! Dissipation only controlled by frame share attenuation in poroelastic (see
+! Morency & Tromp, GJI 2008).
+! attenuation is implemented following the memory variable formulation of
+! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
+! vol. 58(1), p. 110-120 (1993). More details can be found in
+! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a
+! linear
+! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611
+! (1988).
+
+! else
+
+! no attenuation
+ sigma_xx = sigma_xx + &
+ lambdalplus2mul_G*duxdxl + lambdal_G*duydyl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_yy = sigma_yy + &
+ lambdalplus2mul_G*duydyl + lambdal_G*duxdxl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_zz = sigma_zz + &
+ lambdalplus2mul_G*duzdzl + lambdal_G*duxdxl_plus_duydyl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ sigma_xy = sigma_xy + mul_G*duxdyl_plus_duydxl
+ sigma_xz = sigma_xz + mul_G*duzdxl_plus_duxdzl
+ sigma_yz = sigma_yz + mul_G*duzdyl_plus_duydzl
+
+ sigmap = C_biot*duxdxl_plus_duydyl_plus_duzdzl + M_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+
+ ! gets associated normal on GLL point
+ ! (note convention: pointing outwards of poroelastic element)
+ nx = coupling_el_po_normal(1,igll,iface)
+ ny = coupling_el_po_normal(2,igll,iface)
+ nz = coupling_el_po_normal(3,igll,iface)
+
+ ! gets associated, weighted 2D jacobian
+ ! (note: should be the same for poroelastic and elastic element)
+ jacobianw = coupling_el_po_jacobian2Dw(igll,iface)
+
+ ! continuity of displacement and traction on global point
+ !
+ ! note: continuity of displacement is enforced after the velocity update
+! contribution to the solid phase
+ accels_poroelastic(1,iglob_po) = accels_poroelastic(1,iglob_po) + jacobianw*&
+ ( (sigma_xx*nx + sigma_xy*ny + sigma_xz*nz)/2.d0 - phil/tortl*sigmap*nx )
+
+ accels_poroelastic(2,iglob_po) = accels_poroelastic(2,iglob_po) + jacobianw*&
+ ( (sigma_xy*nx + sigma_yy*ny + sigma_yz*nz)/2.d0 - phil/tortl*sigmap*ny )
+
+ accels_poroelastic(3,iglob_po) = accels_poroelastic(3,iglob_po) + jacobianw*&
+ ( (sigma_xz*nx + sigma_yz*ny + sigma_zz*nz)/2.d0 - phil/tortl*sigmap*nz )
+! contribution to the fluid phase
+! w = 0
+
+ enddo ! igll
+
+ endif
+
+ enddo ! iface
+
+end subroutine compute_coupling_poroelastic_el
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -206,9 +206,24 @@
endif
! poroelastic coupling
-! not implemented yet
- !if(POROELASTIC_SIMULATION ) &
- ! call compute_coupling_acoustic_poro()
+ if(POROELASTIC_SIMULATION ) then
+ if( num_coupling_ac_po_faces > 0 ) then
+ if( SIMULATION_TYPE == 1 ) then
+ call compute_coupling_acoustic_po(NSPEC_AB,NGLOB_AB, &
+ ibool,displs_poroelastic,displw_poroelastic, &
+ potential_dot_dot_acoustic, &
+ num_coupling_ac_po_faces, &
+ coupling_ac_po_ispec,coupling_ac_po_ijk, &
+ coupling_ac_po_normal, &
+ coupling_ac_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ else
+ stop 'not implemented yet'
+ endif
+ if( SIMULATION_TYPE == 3 ) &
+ stop 'not implemented yet'
+ endif
+ endif
! sources
call compute_add_sources_acoustic(NSPEC_AB,NGLOB_AB,potential_dot_dot_acoustic, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -167,9 +167,28 @@
! poroelastic coupling
-! not implemented yet
-! if( POROELASTIC_SIMULATION ) &
-! call compute_coupling_elastic_poro()
+ !print *,'entering poro counpling'
+ if( POROELASTIC_SIMULATION ) &
+ call compute_coupling_elastic_po(NSPEC_AB,NGLOB_AB,ibool,&
+ displs_poroelastic,accels_poroelastic,displw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,&
+ displ,accel,kappastore, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store, &
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_coupling_el_po_faces, &
+ coupling_el_po_ispec,coupling_po_el_ispec, &
+ coupling_el_po_ijk,coupling_po_el_ijk, &
+ coupling_el_po_normal, &
+ coupling_el_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ !print *,'ok poro counpling'
! adds source term (single-force/moment-tensor solution)
call compute_add_sources_elastic( NSPEC_AB,NGLOB_AB,accel, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -112,10 +112,28 @@
endif
! elastic coupling
-!chris: TO DO
- if( ELASTIC_SIMULATION ) &
- stop 'elastic-poroelastic simulation not implemented yet'
-! call compute_coupling_poroelastic_el()
+ !print *,'entering elastic counpling'
+ if( ELASTIC_SIMULATION ) then
+ call compute_coupling_poroelastic_el(NSPEC_AB,NGLOB_AB,ibool,&
+ displs_poroelastic,accels_poroelastic,displw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,&
+ displ,accel,kappastore, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store, &
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_coupling_el_po_faces, &
+ coupling_el_po_ispec,coupling_po_el_ispec,&
+ coupling_el_po_ijk,coupling_po_el_ijk, &
+ coupling_el_po_normal, &
+ coupling_el_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+ !print *,'ok elastic counpling'
! adjoint simulations
!chris: TO DO
@@ -127,7 +145,7 @@
! coupling_ac_el_normal, &
! coupling_ac_el_jacobian2Dw, &
! ispec_is_inner,phase_is_inner)
-! endif
+ endif
! adds source term (single-force/moment-tensor solution)
call compute_add_sources_poroelastic( NSPEC_AB,NGLOB_AB, &
@@ -256,5 +274,134 @@
! if (SIMULATION_TYPE == 3) b_velocw_poroelastic(:,:) = b_velocw_poroelastic(:,:) + &
! b_deltatover2*b_accelw_poroelastic(:,:)
+
+
+! elastic coupling
+ !print*, 'entering assembling displacements'
+ if( ELASTIC_SIMULATION ) &
+ call compute_continuity_disp_po_el(NSPEC_AB,NGLOB_AB,ibool,&
+ accel,veloc,&
+ accels_poroelastic,velocs_poroelastic,&
+ accelw_poroelastic,velocw_poroelastic,&
+ rmass,rmass_solid_poroelastic,rmass_fluid_poroelastic,&
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT,&
+ num_coupling_el_po_faces,&
+ coupling_el_po_ispec,coupling_el_po_ijk,&
+ deltatover2)
+ !print*, 'ok assembling displacements'
+
+
end subroutine compute_forces_poroelastic
+
+!-----------------------------------------------------------------------------------------------------------------
+
+subroutine compute_continuity_disp_po_el(NSPEC_AB,NGLOB_AB,ibool,&
+ accel,veloc,&
+ accels_poroelastic,velocs_poroelastic,&
+ accelw_poroelastic,velocw_poroelastic,&
+ rmass,rmass_solid_poroelastic,rmass_fluid_poroelastic,&
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT,&
+ num_coupling_el_po_faces,&
+ coupling_el_po_ispec,coupling_el_po_ijk,&
+ deltatover2)
+!*******************************************************************************
+! assembling the displacements on the elastic-poro boundaries
+!*******************************************************************************
+
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB),intent(inout) :: veloc,accel
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB),intent(inout) :: velocs_poroelastic,accels_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB),intent(inout) :: velocw_poroelastic,accelw_poroelastic
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB),intent(in) :: rmass,rmass_solid_poroelastic,rmass_fluid_poroelastic
+
+! global indexing
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool
+
+ integer :: SIMULATION_TYPE
+ integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+
+! elastic-poroelastic coupling surface
+ integer :: num_coupling_el_po_faces
+ integer :: coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
+ integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
+
+ real(kind=CUSTOM_REAL),intent(in) :: deltatover2
+
+! local variables
+ integer, dimension(NGLOB_AB) :: icount
+ integer :: ispec,i,j,k,iglob,igll,iface
+
+ icount(:)=ZERO
+
+! loops on all coupling faces
+ do iface = 1,num_coupling_el_po_faces
+
+ ! gets corresponding poroelastic spectral element
+ ispec = coupling_el_po_ispec(iface)
+
+ ! loops over common GLL points
+ do igll = 1, NGLLSQUARE
+ i = coupling_el_po_ijk(1,igll,iface)
+ j = coupling_el_po_ijk(2,igll,iface)
+ k = coupling_el_po_ijk(3,igll,iface)
+
+ ! gets global index of this common GLL point
+ iglob = ibool(i,j,k,ispec)
+ icount(iglob) = icount(iglob) + 1
+
+ if(icount(iglob) ==1)then
+
+! recovering original velocities and accelerations on boundaries (elastic side)
+ veloc(1,iglob) = veloc(1,iglob) - deltatover2*accel(1,iglob)
+ veloc(2,iglob) = veloc(2,iglob) - deltatover2*accel(2,iglob)
+ veloc(3,iglob) = veloc(3,iglob) - deltatover2*accel(3,iglob)
+ accel(1,iglob) = accel(1,iglob) / rmass(iglob)
+ accel(2,iglob) = accel(2,iglob) / rmass(iglob)
+ accel(3,iglob) = accel(3,iglob) / rmass(iglob)
+! recovering original velocities and accelerations on boundaries (poro side)
+ velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) - deltatover2*accels_poroelastic(1,iglob)
+ velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) - deltatover2*accels_poroelastic(2,iglob)
+ velocs_poroelastic(3,iglob) = velocs_poroelastic(3,iglob) - deltatover2*accels_poroelastic(3,iglob)
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_solid_poroelastic(iglob)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_solid_poroelastic(iglob)
+ accels_poroelastic(3,iglob) = accels_poroelastic(3,iglob) / rmass_solid_poroelastic(iglob)
+! assembling accelerations
+ accel(1,iglob) = ( accel(1,iglob) + accels_poroelastic(1,iglob) ) / &
+ ( 1.0/rmass(iglob) +1.0/rmass_solid_poroelastic(iglob) )
+ accel(2,iglob) = ( accel(2,iglob) + accels_poroelastic(2,iglob) ) / &
+ ( 1.0/rmass(iglob) +1.0/rmass_solid_poroelastic(iglob) )
+ accel(3,iglob) = ( accel(3,iglob) + accels_poroelastic(3,iglob) ) / &
+ ( 1.0/rmass(iglob) +1.0/rmass_solid_poroelastic(iglob) )
+ accels_poroelastic(1,iglob) = accel(1,iglob)
+ accels_poroelastic(2,iglob) = accel(2,iglob)
+ accels_poroelastic(3,iglob) = accel(3,iglob)
+! updating velocities
+ velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) + deltatover2*accels_poroelastic(1,iglob)
+ velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) + deltatover2*accels_poroelastic(2,iglob)
+ velocs_poroelastic(3,iglob) = velocs_poroelastic(3,iglob) + deltatover2*accels_poroelastic(3,iglob)
+ veloc(1,iglob) = veloc(1,iglob) + deltatover2*accel(1,iglob)
+ veloc(2,iglob) = veloc(2,iglob) + deltatover2*accel(2,iglob)
+ veloc(3,iglob) = veloc(3,iglob) + deltatover2*accel(3,iglob)
+! zeros w
+ accelw_poroelastic(1,iglob) = 0.d0
+ accelw_poroelastic(2,iglob) = 0.d0
+ accelw_poroelastic(3,iglob) = 0.d0
+ velocw_poroelastic(1,iglob) = 0.d0
+ velocw_poroelastic(2,iglob) = 0.d0
+ velocw_poroelastic(3,iglob) = 0.d0
+
+ if(SIMULATION_TYPE == 3) then
+! to do
+ endif
+
+ endif !if(icount(iglob) ==1)
+ enddo ! igll
+ enddo ! iface
+
+end subroutine compute_continuity_disp_po_el
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -78,6 +78,7 @@
if( ACOUSTIC_SIMULATION ) call compute_forces_acoustic()
! elastic solver
+ ! (needs to be done first, before poroelastic one)
if( ELASTIC_SIMULATION ) call compute_forces_elastic()
! poroelastic solver
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -307,11 +307,15 @@
allocate(coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces), &
coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces), &
coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces), &
- coupling_el_po_ispec(num_coupling_el_po_faces),stat=ier)
+ coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces), &
+ coupling_el_po_ispec(num_coupling_el_po_faces), &
+ coupling_po_el_ispec(num_coupling_el_po_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array coupling_el_po_normal etc.'
if( num_coupling_el_po_faces > 0 ) then
read(27) coupling_el_po_ispec
+ read(27) coupling_po_el_ispec
read(27) coupling_el_po_ijk
+ read(27) coupling_po_el_ijk
read(27) coupling_el_po_jacobian2Dw
read(27) coupling_el_po_normal
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -447,8 +447,8 @@
! elastic-poroelastic coupling surface
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: coupling_el_po_normal
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: coupling_el_po_jacobian2Dw
- integer, dimension(:,:,:), allocatable :: coupling_el_po_ijk
- integer, dimension(:), allocatable :: coupling_el_po_ispec
+ integer, dimension(:,:,:), allocatable :: coupling_el_po_ijk,coupling_po_el_ijk
+ integer, dimension(:), allocatable :: coupling_el_po_ispec,coupling_po_el_ispec
integer :: num_coupling_el_po_faces
! material flag
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-11-26 19:01:05 UTC (rev 19244)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-11-29 04:02:32 UTC (rev 19245)
@@ -1022,6 +1022,7 @@
use specfem_par
use specfem_par_elastic
+ use specfem_par_poroelastic
use specfem_par_acoustic
use specfem_par_movie
implicit none
@@ -1104,6 +1105,29 @@
deallocate(div_glob,curl_glob,valency)
+ endif ! elastic
+
+ ! saves full snapshot data to local disk
+ if( POROELASTIC_SIMULATION ) then
+ ! allocate array for single elements
+ allocate(div_glob(NGLOB_AB), &
+ curl_glob(NGLOB_AB), &
+ valency(NGLOB_AB), stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie div and curl'
+ ! calculates divergence and curl of velocity field
+ call wmo_movie_div_curl(NSPEC_AB,NGLOB_AB,velocs_poroelastic, &
+ div_glob,curl_glob,valency, &
+ div,curl_x,curl_y,curl_z, &
+ velocity_x,velocity_y,velocity_z, &
+ ibool,ispec_is_poroelastic, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+ deallocate(div_glob,curl_glob,valency)
+ endif ! poroelastic
+
+
+ if( ELASTIC_SIMULATION .or. POROELASTIC_SIMULATION ) then
+
! writes our divergence
write(outputname,"('/proc',i6.6,'_div_it',i6.6,'.bin')") myrank,it
open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
@@ -1133,9 +1157,10 @@
!write(27) veloc
!close(27)
- endif ! elastic
+ endif
- if( ACOUSTIC_SIMULATION .or. ELASTIC_SIMULATION ) then
+
+ if( ACOUSTIC_SIMULATION .or. ELASTIC_SIMULATION .or. POROELASTIC_SIMULATION) then
write(outputname,"('/proc',i6.6,'_velocity_',a1,'_it',i6.6,'.bin')") myrank,compx,it
open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
if( ier /= 0 ) stop 'error opening file movie output velocity x'
@@ -1170,7 +1195,7 @@
div_glob,curl_glob,valency, &
div,curl_x,curl_y,curl_z, &
velocity_x,velocity_y,velocity_z, &
- ibool,ispec_is_elastic, &
+ ibool,ispec_is, &
hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
@@ -1192,7 +1217,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: div, curl_x, curl_y, curl_z
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: velocity_x,velocity_y,velocity_z
integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB):: ibool
- logical,dimension(NSPEC_AB) :: ispec_is_elastic
+ logical,dimension(NSPEC_AB) :: ispec_is
! array with derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
@@ -1219,7 +1244,7 @@
! loops over elements
do ispec=1,NSPEC_AB
- if( .not. ispec_is_elastic(ispec) ) cycle
+ if( .not. ispec_is(ispec) ) cycle
! calculates divergence and curl of velocity field
do k=1,NGLLZ
More information about the CIG-COMMITS
mailing list