[cig-commits] r19025 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh_SCOTCH generate_databases shared specfem3D
cmorency at geodynamics.org
cmorency at geodynamics.org
Wed Oct 5 19:32:28 PDT 2011
Author: cmorency
Date: 2011-10-05 19:32:28 -0700 (Wed, 05 Oct 2011)
New Revision: 19025
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/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.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/assemble_MPI_vector.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
Log:
Merging of 3D Biot into SPECFEM3D with CUBIT support...
Modify files in src/
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-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -106,6 +106,8 @@
integer :: aniso_flag,idomain_id
double precision :: vp,vs,rho,qmu
+! poroelastic parameters read in a new file
+ double precision :: rhos,rhof,phi,tort,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,eta,mufr
contains
@@ -217,6 +219,9 @@
! vs : S-velocity
! Q_mu : 0=no attenuation
! anisotropy_flag : 0=no anisotropy/ 1,2,.. check with implementation in aniso_model.f90
+ ! Note that when poroelastic material, this file is a dummy except for material_domain_id & material_id,
+ ! and that poroelastic materials are actually read from nummaterial_poroelastic_file, because CUBIT
+ ! cannot support more than 10 attributes
count_def_mat = 0
count_undef_mat = 0
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file',&
@@ -247,7 +252,7 @@
print*,' bigger than defined materials in nummaterial_velocity_file:',count_def_mat
stop 'error materials'
endif
- allocate(mat_prop(6,count_def_mat),stat=ier)
+ allocate(mat_prop(16,count_def_mat),stat=ier)
if( ier /= 0 ) stop 'error allocating array mat_prop'
allocate(undef_mat_prop(6,count_undef_mat),stat=ier)
if( ier /= 0 ) stop 'error allocating array undef_mat_prop'
@@ -259,6 +264,24 @@
status='old', form='formatted', iostat=ier)
if( ier /= 0 ) stop 'error opening nummaterial_velocity_file'
+ ! modif to read poro parameters, added if loop on idomain_id
+ ! note: format of nummaterial_poroelastic_file located in MESH must be
+ !
+ ! #(1)rhos,#(2)rhof,#(3)phi,#(4)tort,#(5)kxx,#(6)kxy,#(7)kxz,#(8)kyy,#(9)kyz,#(10)kzz,
+ ! #(11)kappas,#(12)kappaf,#(13)kappafr,#(14)eta,#(15)mufr
+ !
+ ! where
+ ! rhos, rhof : solid & fluid density
+ ! phi : porosity
+ ! tort : tortuosity
+ ! k : permeability tensor
+ ! kappas, kappaf, kappafr : solid, fluid and frame bulk moduli
+ ! eta : fluid viscosity
+ ! mufr : frame shear modulus
+ open(unit=97, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_poroelastic_file', &
+ status='old', form='formatted', iostat=ier)
+ if( ier /= 0 ) stop 'error opening nummaterial_poroelastic_file'
+
! note: entries in nummaterial_velocity_file can be an unsorted list of all
! defined materials (material_id > 0) and undefined materials (material_id < 0 )
do imat=1,count_def_mat
@@ -282,6 +305,8 @@
! checks material_id bounds
if(num_mat < 1 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
+ if(idomain_id == 1 .or. idomain_id == 2) then ! material is elastic or acoustic
+
!read(98,*) num_mat, mat_prop(1,num_mat),mat_prop(2,num_mat),&
! mat_prop(3,num_mat),mat_prop(4,num_mat),mat_prop(5,num_mat)
mat_prop(1,num_mat) = rho
@@ -291,6 +316,28 @@
mat_prop(5,num_mat) = aniso_flag
mat_prop(6,num_mat) = idomain_id
+ else ! material is poroelastic
+
+ read(97,*) rhos,rhof,phi,tort,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,eta,mufr
+ mat_prop(1,num_mat) = rhos
+ mat_prop(2,num_mat) = rhof
+ mat_prop(3,num_mat) = phi
+ mat_prop(4,num_mat) = tort
+ mat_prop(5,num_mat) = eta
+ mat_prop(6,num_mat) = idomain_id
+ mat_prop(7,num_mat) = kxx
+ mat_prop(8,num_mat) = kxy
+ mat_prop(9,num_mat) = kxz
+ mat_prop(10,num_mat) = kyy
+ mat_prop(11,num_mat) = kyz
+ mat_prop(12,num_mat) = kzz
+ mat_prop(13,num_mat) = kappas
+ mat_prop(14,num_mat) = kappaf
+ mat_prop(15,num_mat) = kappafr
+ mat_prop(16,num_mat) = mufr
+
+ endif !if(idomain_id == 1 .or. idomain_id == 2)
+
end do
! reads in undefined material properties
@@ -370,6 +417,7 @@
endif
endif
end do
+ close(97)
close(98)
@@ -638,7 +686,7 @@
num_material(:) = mat(1,:)
! in case of acoustic/elastic simulation, weights elements accordingly
- call acoustic_elastic_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
+ call acoustic_elastic_poroelastic_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
num_material,mat_prop,undef_mat_prop)
deallocate(num_material)
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-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -40,10 +40,12 @@
! very large and very small values
double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
-! acoustic-elastic load balancing:
+! acoustic-elastic-poroelastic load balancing:
! assumes that elastic at least ~4 times more expensive than acoustic
+! assumes that poroelastic at least ~8 times more expensive than acoustic
integer, parameter :: ACOUSTIC_LOAD = 1
integer, parameter :: ELASTIC_LOAD = 4
+ integer, parameter :: POROELASTIC_LOAD = 8
! include './constants_decompose_mesh_SCOTCH.h'
@@ -282,7 +284,7 @@
! 1/ first element, 2/ second element, 3/ number of common nodes, 4/ first node,
! 5/ second node, if relevant.
- ! interface ignores acoustic and elastic elements
+ ! interface ignores acoustic, elastic and poroelastic elements
! Elements with undefined material are considered as elastic elements.
!--------------------------------------------------
@@ -607,7 +609,7 @@
integer, intent(in) :: IIN_database
integer, intent(in) :: count_def_mat,count_undef_mat
- double precision, dimension(6,count_def_mat) :: mat_prop
+ double precision, dimension(16,count_def_mat) :: mat_prop
character (len=30), dimension(6,count_undef_mat) :: undef_mat_prop
integer :: i
@@ -615,12 +617,17 @@
do i = 1, count_def_mat
! database material definition
!
- ! format: #rho #vp #vs #Q_flag #anisotropy_flag #domain_id
+ ! format: #rho #vp #vs #Q_flag #anisotropy_flag #domain_id for elastic/acoustic
+ ! format: #rhos,#rhof,#phi,#tort,#kxx,#domain_id,#kxy,#kxz,#kyy,#kyz,#kzz,
+ ! #kappas,#kappaf,#kappafr,#eta,#mufr for poroelastic
!
! (note that this order of the properties is different than the input in nummaterial_velocity_file)
!
write(IIN_database,*) mat_prop(1,i), mat_prop(2,i), mat_prop(3,i), &
- mat_prop(4,i), mat_prop(5,i), mat_prop(6,i)
+ mat_prop(4,i), mat_prop(5,i), mat_prop(6,i), &
+ mat_prop(7,i), mat_prop(8,i), mat_prop(9,i), &
+ mat_prop(10,i), mat_prop(11,i), mat_prop(12,i), &
+ mat_prop(13,i), mat_prop(14,i), mat_prop(15,i), mat_prop(16,i)
end do
do i = 1, count_undef_mat
write(IIN_database,*) trim(undef_mat_prop(1,i)),' ',trim(undef_mat_prop(2,i)),' ', &
@@ -1230,16 +1237,17 @@
!--------------------------------------------------
- ! loading : sets weights for acoustic/elastic elements to account for different
+ ! loading : sets weights for acoustic/elastic/poroelastic elements to account for different
! expensive calculations in specfem simulations
!--------------------------------------------------
- subroutine acoustic_elastic_load (elmnts_load,nelmnts,count_def_mat,count_undef_mat, &
+ subroutine acoustic_elastic_poroelastic_load (elmnts_load,nelmnts,count_def_mat,count_undef_mat, &
num_material,mat_prop,undef_mat_prop)
!
! note:
! acoustic material = domainID 1 (stored in mat_prop(6,..) )
! elastic material = domainID 2
+ ! poroelastic material = domainID 3
!
implicit none
@@ -1251,18 +1259,19 @@
! materials
integer, dimension(1:nelmnts), intent(in) :: num_material
- double precision, dimension(6,count_def_mat),intent(in) :: mat_prop
+ double precision, dimension(16,count_def_mat),intent(in) :: mat_prop
character (len=30), dimension(6,count_undef_mat),intent(in) :: undef_mat_prop
! local parameters
- logical, dimension(-count_undef_mat:count_def_mat) :: is_acoustic, is_elastic
+ logical, dimension(-count_undef_mat:count_def_mat) :: is_acoustic, is_elastic, is_poroelastic
integer :: i,el,idomain_id
! initializes flags
is_acoustic(:) = .false.
is_elastic(:) = .false.
+ is_poroelastic(:) = .false.
- ! sets acoustic/elastic flags for defined materials
+ ! sets acoustic/elastic/poroelastic flags for defined materials
do i = 1, count_def_mat
idomain_id = mat_prop(6,i)
! acoustic material has idomain_id 1
@@ -1273,6 +1282,10 @@
if (idomain_id == 2 ) then
is_elastic(i) = .true.
endif
+ ! poroelastic material has idomain_id 3
+ if (idomain_id == 3 ) then
+ is_poroelastic(i) = .true.
+ endif
enddo
! sets acoustic/elastic flags for undefined materials
@@ -1299,9 +1312,13 @@
if ( is_elastic(num_material(el+1)) ) then
elmnts_load(el+1) = elmnts_load(el+1)*ELASTIC_LOAD
endif
+ ! poroelastic element (very expensive)
+ if ( is_poroelastic(num_material(el+1)) ) then
+ elmnts_load(el+1) = elmnts_load(el+1)*POROELASTIC_LOAD
+ endif
enddo
- end subroutine acoustic_elastic_load
+ end subroutine acoustic_elastic_poroelastic_load
!--------------------------------------------------
@@ -1321,7 +1338,7 @@
integer, dimension(1:nelmnts), intent(in) :: num_material
- double precision, dimension(6,nb_materials),intent(in) :: mat_prop
+ double precision, dimension(16,nb_materials),intent(in) :: mat_prop
integer, dimension(0:nelmnts-1) :: part
integer, dimension(0:esize*nelmnts-1) :: elmnts
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -42,6 +42,7 @@
double precision :: weight
real(kind=CUSTOM_REAL) :: jacobianl
integer :: ispec,i,j,k,iglob,ier
+ real(kind=CUSTOM_REAL) :: rho_s, rho_f, rho_bar, phi, tort
! allocates memory
allocate(rmass(nglob),stat=ier); if(ier /= 0) stop 'error in allocate rmass'
@@ -90,31 +91,30 @@
! poroelastic mass matrices
if( ispec_is_poroelastic(ispec) ) then
- stop 'poroelastic mass matrices not implemented yet'
+ rho_s = rhoarraystore(1,i,j,k,ispec)
+ rho_f = rhoarraystore(2,i,j,k,ispec)
+ phi = phistore(i,j,k,ispec)
+ tort = tortstore(i,j,k,ispec)
+ rho_bar = (1._CUSTOM_REAL-phi)*rho_s + phi*rho_f
- !rho_solid = density(1,kmato(ispec))
- !rho_fluid = density(2,kmato(ispec))
- !phi = porosity(kmato(ispec))
- !tort = tortuosity(kmato(ispec))
- !rho_bar = (1._CUSTOM_REAL-phil)*rhol_s + phil*rhol_f
- !
- !if(CUSTOM_REAL == SIZE_REAL) then
- ! ! for the solid mass matrix
- ! rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
- ! sngl( dble(jacobianl) * weight * dble(rho_bar - phi*rho_fluid/tort) )
- !
- ! ! for the fluid mass matrix
- ! rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
- ! sngl( dble(jacobianl) * weight * dble(rho_bar*rho_fluid*tort - &
- ! phi*rho_fluid*rho_fluid)/dble(rho_bar*phi) )
- !else
- ! rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
- ! jacobianl * weight * (rho_bar - phi*rho_fluid/tort)
- !
- ! rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
- ! jacobianl * weight * (rho_bar*rho_fluid*tort - &
- ! phi*rho_fluid*rho_fluid) / (rho_bar*phi)
- !endif
+ if(CUSTOM_REAL == SIZE_REAL) then
+ ! for the solid mass matrix
+ rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(rho_bar - phi*rho_f/tort) )
+
+ ! for the fluid mass matrix
+ rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(rho_bar*rho_f*tort - &
+ phi*rho_f*rho_f)/dble(rho_bar*phi) )
+ else
+ rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
+ jacobianl * weight * (rho_bar - phi*rho_f/tort)
+
+ rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
+ jacobianl * weight * (rho_bar*rho_f*tort - &
+ phi*rho_f*rho_f) / (rho_bar*phi)
+ endif
+
endif
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -49,6 +49,11 @@
! for model density, kappa, mu
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
+! for poroelastic model
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:), allocatable :: etastore,phistore,tortstore
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:,:), allocatable :: rhoarraystore,kappaarraystore,permstore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vpI,rho_vpII,rho_vsI
+
! mass matrix
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass,rmass_acoustic,&
rmass_solid_poroelastic,rmass_fluid_poroelastic
@@ -86,6 +91,20 @@
integer, dimension(:), allocatable :: coupling_ac_el_ispec
integer :: num_coupling_ac_el_faces
+! acoustic-poroelastic coupling surface
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: coupling_ac_po_normal
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: coupling_ac_po_jacobian2Dw
+ integer, dimension(:,:,:), allocatable :: coupling_ac_po_ijk
+ integer, dimension(:), allocatable :: coupling_ac_po_ispec
+ integer :: num_coupling_ac_po_faces
+
+! 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 :: num_coupling_el_po_faces
+
! for stacey
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
@@ -166,7 +185,7 @@
! material properties
integer :: nmat_ext_mesh,nundefMat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
+ double precision, dimension(16,nmat_ext_mesh) :: materials_ext_mesh
character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
! double precision, external :: materials_ext_mesh
@@ -307,10 +326,10 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
nspec2D_bottom,nspec2D_top)
-! sets up acoustic-elastic coupling surfaces
+! sets up acoustic-elastic-poroelastic coupling surfaces
call sync_all()
if( myrank == 0) then
- write(IMAIN,*) ' ...detecting acoustic-elastic surfaces '
+ write(IMAIN,*) ' ...detecting acoustic-elastic-poroelastic surfaces '
endif
call get_coupling_surfaces(myrank, &
nspec,nglob,ibool,NPROC, &
@@ -346,6 +365,8 @@
gammaxstore,gammaystore,gammazstore, &
jacobianstore, rho_vp,rho_vs,qmu_attenuation_store, &
rhostore,kappastore,mustore, &
+ rhoarraystore,kappaarraystore,etastore,phistore,tortstore,permstore, &
+ rho_vpI,rho_vpII,rho_vsI, &
rmass,rmass_acoustic,rmass_solid_poroelastic,rmass_fluid_poroelastic, &
OCEANS,rmass_ocean_load,NGLOB_OCEAN,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
abs_boundary_normal,abs_boundary_jacobian2Dw, &
@@ -354,6 +375,10 @@
free_surface_ijk,free_surface_ispec,num_free_surface_faces, &
coupling_ac_el_normal,coupling_ac_el_jacobian2Dw, &
coupling_ac_el_ijk,coupling_ac_el_ispec,num_coupling_ac_el_faces, &
+ coupling_ac_po_normal,coupling_ac_po_jacobian2Dw, &
+ 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, &
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,ANISOTROPY,NSPEC_ANISO, &
@@ -370,10 +395,12 @@
! checks the mesh, stability and resolved period
call sync_all()
+!chris: check for poro: At the moment cpI & cpII are for eta=0
call check_mesh_resolution(myrank,nspec,nglob,ibool,&
xstore_dummy,ystore_dummy,zstore_dummy, &
kappastore,mustore,rho_vp,rho_vs, &
- -1.0d0, model_speed_max,min_resolved_period )
+ -1.0d0, model_speed_max,min_resolved_period, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI )
! saves binary mesh files for attenuation
if( ATTENUATION ) then
@@ -414,6 +441,8 @@
gammaxstore,gammaystore,gammazstore)
deallocate(jacobianstore,qmu_attenuation_store)
deallocate(kappastore,mustore,rho_vp,rho_vs)
+ deallocate(rho_vpI,rho_vpII,rho_vsI)
+ deallocate(rhoarraystore,kappaarraystore,etastore,phistore,tortstore,permstore)
end subroutine create_regions_mesh_ext
@@ -507,6 +536,18 @@
!vsstore(NGLLX,NGLLY,NGLLZ,nspec),
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+! array with poroelastic model
+ allocate(rhoarraystore(2,NGLLX,NGLLY,NGLLZ,nspec), &
+ kappaarraystore(3,NGLLX,NGLLY,NGLLZ,nspec), &
+ etastore(NGLLX,NGLLY,NGLLZ,nspec), &
+ tortstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ phistore(NGLLX,NGLLY,NGLLZ,nspec), &
+ rho_vpI(NGLLX,NGLLY,NGLLZ,nspec), &
+ rho_vpII(NGLLX,NGLLY,NGLLZ,nspec), &
+ rho_vsI(NGLLX,NGLLY,NGLLZ,nspec), &
+ permstore(6,NGLLX,NGLLY,NGLLZ,nspec), stat=ier)
+ if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+
! arrays with mesh parameters
allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec), &
xiystore(NGLLX,NGLLY,NGLLZ,nspec), &
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -438,8 +438,6 @@
if(myrank == 0) then
-! chris: I am not sure if we should suppress the following. topography should appear in the external mesh
-! leave it for now
write(IMAIN,*)
if(SUPPRESS_UTM_PROJECTION) then
@@ -576,13 +574,20 @@
call sync_all()
! read materials' physical properties
+! added poroelastic properties and filled with 0 the last 10 entries for elastic/acoustic
read(IIN,*) nmat_ext_mesh, nundefMat_ext_mesh
- allocate(materials_ext_mesh(6,nmat_ext_mesh),stat=ier)
+ allocate(materials_ext_mesh(16,nmat_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating array materials_ext_mesh'
do imat = 1, nmat_ext_mesh
- ! format: #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag #(6) material_domain_id
+ ! ielastic/acoustic format: #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag #(6) material_domain_id
+ ! and remaining entries are filled with zeros
+ ! poroelastic format: rhos,rhof,phi,tort,eta,material_domain_id,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,mufr
read(IIN,*) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
- materials_ext_mesh(4,imat), materials_ext_mesh(5,imat), materials_ext_mesh(6,imat)
+ materials_ext_mesh(4,imat), materials_ext_mesh(5,imat), materials_ext_mesh(6,imat), &
+ materials_ext_mesh(7,imat), materials_ext_mesh(8,imat), materials_ext_mesh(9,imat), &
+ materials_ext_mesh(10,imat), materials_ext_mesh(11,imat), materials_ext_mesh(12,imat), &
+ materials_ext_mesh(13,imat), materials_ext_mesh(14,imat), materials_ext_mesh(15,imat), &
+ materials_ext_mesh(16,imat)
! output
!print*,'materials:',materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -29,10 +29,10 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
my_neighbours_ext_mesh)
+
+! determines coupling surface for acoustic-elastic-poroelastic domains
-! determines coupling surface for acoustic-elastic domains
-
- use create_regions_mesh_ext_par
+ use create_regions_mesh_ext_par
implicit none
! number of spectral elements in each block
@@ -50,7 +50,7 @@
! local parameters
! (assumes NGLLX=NGLLY=NGLLZ)
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
real(kind=CUSTOM_REAL),dimension(:,:,:),allocatable :: tmp_normal
@@ -61,15 +61,16 @@
integer,dimension(NGNOD2D) :: iglob_corners_ref !,iglob_corners
integer :: ispec,i,j,k,igll,ier,iglob
- integer :: inum,iface_ref,icorner,iglob_midpoint ! iface,ispec_neighbor
- integer :: count_elastic,count_acoustic
-
+ integer :: inum,inum_ac,inum_el,inum_po,iface_ref,icorner,iglob_midpoint ! iface,ispec_neighbor
+ integer :: inum_ac_el,inum_el_po,inum_ac_po
+ integer :: count_elastic,count_acoustic,count_poroelastic
+
! mpi interface communication
- integer, dimension(:), allocatable :: elastic_flag,acoustic_flag,test_flag
+ integer, dimension(:), allocatable :: elastic_flag,acoustic_flag,poroelastic_flag,test_flag
integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
integer :: max_nibool_interfaces_ext_mesh
- logical, dimension(:), allocatable :: mask_ibool
-
+ logical, dimension(:), allocatable :: mask_ibool_ac_el,mask_ibool_ac_po,mask_ibool_el_po
+
! corners indices of reference cube faces
integer,dimension(3,4),parameter :: iface1_corner_ijk = &
reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
@@ -82,21 +83,21 @@
integer,dimension(3,4),parameter :: iface5_corner_ijk = &
reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
integer,dimension(3,4),parameter :: iface6_corner_ijk = &
- reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
+ reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
iface3_corner_ijk,iface4_corner_ijk, &
iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
+ ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
integer,dimension(3,6),parameter :: iface_all_midpointijk = &
- reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
+ reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
-
+
! test vtk output
!integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
!character(len=256):: prname_file
-
-! allocates temporary arrays
+
+! allocates temporary arrays
allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6),stat=ier)
if( ier /= 0 ) stop 'error allocating array tmp_normal'
allocate(tmp_jacobian2Dw(NGLLSQUARE,nspec*6),stat=ier)
@@ -109,49 +110,62 @@
tmp_ijk(:,:,:) = 0
tmp_normal(:,:,:) = 0.0
tmp_jacobian2Dw(:,:) = 0.0
-
- ! sets flags for acoustic / elastic on global points
+
+ ! sets flags for acoustic / elastic /poroelastic on global points
allocate(elastic_flag(nglob),stat=ier)
if( ier /= 0 ) stop 'error allocating array elastic_flag'
allocate(acoustic_flag(nglob),stat=ier)
if( ier /= 0 ) stop 'error allocating array acoustic_flag'
+ allocate(poroelastic_flag(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array poroelastic_flag'
allocate(test_flag(nglob),stat=ier)
if( ier /= 0 ) stop 'error allocating array test_flag'
- allocate(mask_ibool(nglob),stat=ier)
+ allocate(mask_ibool_ac_el(nglob),stat=ier)
+ allocate(mask_ibool_ac_po(nglob),stat=ier)
+ allocate(mask_ibool_el_po(nglob),stat=ier)
if( ier /= 0 ) stop 'error allocating array mask_ibool'
elastic_flag(:) = 0
acoustic_flag(:) = 0
+ poroelastic_flag(:) = 0
test_flag(:) = 0
count_elastic = 0
count_acoustic = 0
+ count_poroelastic = 0
do ispec = 1, nspec
! counts elements
if( ispec_is_elastic(ispec) ) count_elastic = count_elastic + 1
if( ispec_is_acoustic(ispec) ) count_acoustic = count_acoustic + 1
-
+ if( ispec_is_poroelastic(ispec) ) count_poroelastic = count_poroelastic + 1
+
! sets flags on global points
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
! global index
- iglob = ibool(i,j,k,ispec)
+ iglob = ibool(i,j,k,ispec)
! sets elastic flag
if( ispec_is_elastic(ispec) ) elastic_flag(iglob) = myrank+1
! sets acoustic flag
if( ispec_is_acoustic(ispec) ) acoustic_flag(iglob) = myrank+1
+ ! sets poroelastic flag
+ if( ispec_is_poroelastic(ispec) ) poroelastic_flag(iglob) = myrank+1
! sets test flag
test_flag(iglob) = myrank+1
enddo
enddo
enddo
enddo
- call sum_all_i(count_acoustic,inum)
+ call sum_all_i(count_acoustic,inum_ac)
if( myrank == 0 ) then
- write(IMAIN,*) ' total acoustic elements:',inum
- endif
- call sum_all_i(count_elastic,inum)
+ write(IMAIN,*) ' total acoustic elements:',inum_ac
+ endif
+ call sum_all_i(count_elastic,inum_el)
if( myrank == 0 ) then
- write(IMAIN,*) ' total elastic elements :',inum
+ write(IMAIN,*) ' total elastic elements :',inum_el
+ endif
+ call sum_all_i(count_poroelastic,inum_po)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' total poroelastic elements :',inum_po
endif
! collects contributions from different MPI partitions
@@ -161,7 +175,7 @@
if( ier /= 0 ) stop 'error allocating array ibool_interfaces_ext_mesh_dummy'
do i = 1, num_interfaces_ext_mesh
ibool_interfaces_ext_mesh_dummy(:,i) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,i)
- enddo
+ enddo
! sums elastic flags
call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,elastic_flag, &
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
@@ -172,21 +186,29 @@
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
my_neighbours_ext_mesh)
-
+ ! sums poroelastic flags
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,poroelastic_flag, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
+ my_neighbours_ext_mesh)
! sums test flags
call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,test_flag, &
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
my_neighbours_ext_mesh)
- ! loops over all element faces and
+!----------------------
+! acoustic-elastic
+!----------------------
+ inum_ac_el = 0
+ !if (inum_el >0 .and. inum_ac >0) then
+ ! loops over all element faces and
! counts number of coupling faces between acoustic and elastic elements
- mask_ibool(:) = .false.
- inum = 0
+ mask_ibool_ac_el(:) = .false.
do ispec=1,nspec
! loops over each face
- do iface_ref= 1, 6
+ do iface_ref= 1, 6
! takes indices of corners of reference face
do icorner = 1,NGNOD2D
@@ -199,15 +221,15 @@
! reference corner coordinates
xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
- zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
+ zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
enddo
-
+
! checks if face has acoustic side
if( acoustic_flag( iglob_corners_ref(1) ) >= 1 .and. &
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 elastic side
+ acoustic_flag( iglob_corners_ref(4) ) >= 1) then
+ ! checks if face is has an elastic side
if( elastic_flag( iglob_corners_ref(1) ) >= 1 .and. &
elastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
elastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
@@ -216,22 +238,21 @@
! 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)
+ 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. elastic_flag(iglob_midpoint) >= 1) ) then
+ if( .not. mask_ibool_ac_el(iglob_midpoint) ) then
! gets face GLL points i,j,k indices from element face
call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
-
+
! takes each element face only once, if it lies on an MPI interface
! note: this is not exactly load balanced
! lowest rank process collects as many faces as possible, second lowest as so forth
if( (test_flag(iglob_midpoint) == myrank+1) .or. &
(test_flag(iglob_midpoint) > 2*(myrank+1)) ) then
-
+
! gets face GLL 2Djacobian, weighted from element face
call get_jacobian_boundary_face(myrank,nspec, &
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
@@ -248,31 +269,31 @@
ibool,nspec,nglob, &
xstore_dummy,ystore_dummy,zstore_dummy, &
normal_face(:,i,j) )
- ! makes sure that it always points away from acoustic element,
+ ! makes sure that it always points away from acoustic element,
! otherwise switch direction
if( ispec_is_elastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
enddo
enddo
! stores informations about this face
- inum = inum + 1
- tmp_ispec(inum) = ispec
+ inum_ac_el = inum_ac_el + 1
+ tmp_ispec(inum_ac_el) = ispec
igll = 0
do j=1,NGLLY
do i=1,NGLLX
! adds all gll points on this face
igll = igll + 1
-
+
! do we need to store local i,j,k,ispec info? or only global indices iglob?
- tmp_ijk(:,igll,inum) = ijk_face(:,i,j)
-
+ tmp_ijk(:,igll,inum_ac_el) = 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)
-
+ tmp_jacobian2Dw(igll,inum_ac_el) = jacobian2Dw_face(i,j)
+ tmp_normal(:,igll,inum_ac_el) = 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.
+ mask_ibool_ac_el(iglob) = .true.
enddo
enddo
else
@@ -280,21 +301,23 @@
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.
+ mask_ibool_ac_el(iglob) = .true.
enddo
enddo
endif ! test_flag
- endif ! mask_ibool
+ endif ! mask_ibool
endif ! elastic_flag
endif ! acoustic_flag
enddo ! iface_ref
enddo ! ispec
+
+ !endif !if (count_elastic >0 .and. count_acoustic >0)
-! stores completed coupling face informations
-!
-! note: no need to store material parameters on these coupling points
+! stores completed coupling face informations
+!
+! note: no need to store material parameters on these coupling points
! for acoustic-elastic interface
- num_coupling_ac_el_faces = inum
+ num_coupling_ac_el_faces = inum_ac_el
allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array coupling_ac_el_normal'
allocate(coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces),stat=ier)
@@ -307,15 +330,252 @@
coupling_ac_el_normal(:,:,inum) = tmp_normal(:,:,inum)
coupling_ac_el_jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
coupling_ac_el_ijk(:,:,inum) = tmp_ijk(:,:,inum)
- coupling_ac_el_ispec(inum) = tmp_ispec(inum)
+ coupling_ac_el_ispec(inum) = tmp_ispec(inum)
enddo
! user output
- call sum_all_i(num_coupling_ac_el_faces,inum)
+! makes sure processes are synchronized
+ call sum_all_i(num_coupling_ac_el_faces,inum_ac_el)
if( myrank == 0 ) then
write(IMAIN,*) ' acoustic-elastic coupling:'
- write(IMAIN,*) ' total number of faces = ',inum
+ write(IMAIN,*) ' total number of faces = ',inum_ac_el
+ endif
+
+
+!----------------------
+! acoustic-poroelastic
+!----------------------
+ tmp_ispec(:) = 0
+ tmp_ijk(:,:,:) = 0
+ tmp_normal(:,:,:) = 0.0
+ tmp_jacobian2Dw(:,:) = 0.0
+ inum_ac_po = 0
+ !if (inum_po >0 .and. inum_ac >0) then
+ ! loops over all element faces and
+ ! counts number of coupling faces between acoustic and poroelastic elements
+ do ispec=1,nspec
+
+ if(ispec_is_poroelastic(ispec)) then
+
+ ! loops over each face
+ do iface_ref= 1, 6
+
+ ! takes indices of corners of reference face
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface_ref)
+ j = iface_all_corner_ijk(2,icorner,iface_ref)
+ k = iface_all_corner_ijk(3,icorner,iface_ref)
+ ! global reference indices
+ iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
+
+ ! reference corner coordinates
+ xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
+ ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
+ zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
+ enddo
+
+ ! checks if face has acoustic side
+ if( acoustic_flag( iglob_corners_ref(1) ) >= 1 .and. &
+ acoustic_flag( iglob_corners_ref(2) ) >= 1 .and. &
+ acoustic_flag( iglob_corners_ref(3) ) >= 1 .and. &
+ acoustic_flag( iglob_corners_ref(4) ) >= 1) then
+
+
+ ! gets face GLL points i,j,k indices from element face
+ call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
+
+ ! gets face GLL 2Djacobian, weighted from element face
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+
+ ! normal convention: points away from acoustic, reference element
+ ! switch normal direction if necessary
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! directs normals such that they point outwards of element
+ call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ ! makes sure that it always points away from acoustic element,
+ ! otherwise switch direction
+ !if( ispec_is_poroelastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
+ normal_face(:,i,j) = - normal_face(:,i,j)
+ enddo
+ enddo
+
+ ! stores informations about this face
+ inum_ac_po = inum_ac_po + 1
+ tmp_ispec(inum_ac_po) = ispec
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds all gll points on this face
+ igll = igll + 1
+
+ ! do we need to store local i,j,k,ispec info? or only global indices iglob?
+ tmp_ijk(:,igll,inum_ac_po) = ijk_face(:,i,j)
+
+ ! stores weighted jacobian and normals
+ tmp_jacobian2Dw(igll,inum_ac_po) = jacobian2Dw_face(i,j)
+ tmp_normal(:,igll,inum_ac_po) = normal_face(:,i,j)
+
+ enddo
+ enddo
+ endif ! acoustic_flag
+ enddo ! iface_ref
+ endif ! ispec_is_poroelastic
+ enddo ! ispec
+
+ !endif !if (count_poroelastic >0 .and. count_acoustic >0)
+
+! stores completed coupling face informations
+!
+! 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_ac_po
+ 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'
+ allocate(coupling_ac_po_jacobian2Dw(NGLLSQUARE,num_coupling_ac_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_po_jacobian2Dw'
+ allocate(coupling_ac_po_ijk(3,NGLLSQUARE,num_coupling_ac_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_po_ijk'
+ allocate(coupling_ac_po_ispec(num_coupling_ac_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_po_ispec'
+ do inum = 1,num_coupling_ac_po_faces
+ coupling_ac_po_normal(:,:,inum) = tmp_normal(:,:,inum)
+ coupling_ac_po_jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
+ coupling_ac_po_ijk(:,:,inum) = tmp_ijk(:,:,inum)
+ coupling_ac_po_ispec(inum) = tmp_ispec(inum)
+ enddo
+
+! user output
+! makes sure processes are synchronized
+ call sum_all_i(num_coupling_ac_po_faces,inum_ac_po)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' acoustic-poroelastic coupling:'
+ write(IMAIN,*) ' total number of faces = ',inum_ac_po
endif
+
+!----------------------
+! elastic-poroelastic
+!----------------------
+ tmp_ispec(:) = 0
+ tmp_ijk(:,:,:) = 0
+ tmp_normal(:,:,:) = 0.0
+ tmp_jacobian2Dw(:,:) = 0.0
+ inum_el_po = 0
+ !if (inum_el >0 .and. inum_po >0) then
+ ! loops over all element faces and
+ ! counts number of coupling faces between elastic and poroelastic elements
+ do ispec=1,nspec
+
+ if(ispec_is_poroelastic(ispec)) then
+
+ ! loops over each face
+ do iface_ref= 1, 6
+
+ ! takes indices of corners of reference face
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface_ref)
+ j = iface_all_corner_ijk(2,icorner,iface_ref)
+ k = iface_all_corner_ijk(3,icorner,iface_ref)
+ ! global reference indices
+ iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
+
+ ! reference corner coordinates
+ xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
+ ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
+ zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
+ enddo
+
+ ! checks if face has elastic side
+ if( elastic_flag( iglob_corners_ref(1) ) >= 1 .and. &
+ elastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
+ elastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
+ elastic_flag( iglob_corners_ref(4) ) >= 1) then
+
+ ! gets face GLL points i,j,k indices from element face
+ call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
+
+ ! gets face GLL 2Djacobian, weighted from element face
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+
+ ! normal convention: points away from poroelastic, reference element
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! 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) )
+ enddo
+ enddo
+
+ ! stores informations about this face
+ inum_el_po = inum_el_po + 1
+ tmp_ispec(inum_el_po) = ispec
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds all gll points on this face
+ igll = igll + 1
+
+ ! do we need to store local i,j,k,ispec info? or only global indices iglob?
+ tmp_ijk(:,igll,inum_el_po) = ijk_face(:,i,j)
+
+ ! stores weighted jacobian and normals
+ tmp_jacobian2Dw(igll,inum_el_po) = jacobian2Dw_face(i,j)
+ tmp_normal(:,igll,inum_el_po) = normal_face(:,i,j)
+
+ enddo
+ enddo
+ endif ! elastic_flag
+ enddo ! iface_ref
+ endif ! ispec_is_poroelastic
+ enddo ! ispec
+
+ !endif !if (count_elastic >0 .and. count_poroelastic >0)
+
+! stores completed coupling face informations
+!
+! 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_el_po
+ 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'
+ allocate(coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces),stat=ier)
+ 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_el_po_ispec(num_coupling_el_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_el_po_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_el_po_ispec(inum) = tmp_ispec(inum)
+ enddo
+
+! user output
+! makes sure processes are synchronized
+ call sum_all_i(num_coupling_el_po_faces,inum_el_po)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' elastic-poroelastic coupling:'
+ write(IMAIN,*) ' total number of faces = ',inum_el_po
+ endif
+
+
+
end subroutine get_coupling_surfaces
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -42,7 +42,7 @@
integer :: nmat_ext_mesh,nundefMat_ext_mesh
integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
+ double precision, dimension(16,nmat_ext_mesh) :: materials_ext_mesh
character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
! anisotropy
@@ -52,6 +52,10 @@
real(kind=CUSTOM_REAL) :: vp,vs,rho,qmu_atten
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
+ real(kind=CUSTOM_REAL) :: kappa_s,kappa_f,kappa_fr,mu_fr,rho_s,rho_f,phi,tort,eta_f, &
+ kxx,kxy,kxz,kyy,kyz,kzz,rho_bar
+ real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,H_biot,M_biot,C_biot,D_biot, &
+ afactor,bfactor,cfactor
integer :: ispec,i,j,k,iundef,ier
integer :: iflag,flag_below,flag_above
integer :: iflag_aniso,idomain_id,imaterial_id
@@ -103,8 +107,14 @@
! check if the material is known or unknown
if( imaterial_id > 0) then
- ! gets velocity model as specified by (cubit) mesh files
+ ! gets velocity model as specified by (cubit) mesh files for elastic & acoustic
+ ! or from nummaterial_poroelastic_file for poroelastic (too many arguments for cubit)
+ ! material domain_id
+ idomain_id = materials_ext_mesh(6,imaterial_id)
+
+ if(idomain_id == IDOMAIN_ACOUSTIC .or. idomain_id == IDOMAIN_ELASTIC) then ! elastic or acoustic
+
! density
! materials_ext_mesh format:
! #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
@@ -120,9 +130,30 @@
! anisotropy
iflag_aniso = materials_ext_mesh(5,imaterial_id)
- ! material domain_id
- idomain_id = materials_ext_mesh(6,imaterial_id)
+ else ! poroelastic
+ ! materials_ext_mesh format:
+ ! rhos,rhof,phi,tort,eta,domain_id,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,mufr
+ ! solid properties
+ rho_s = materials_ext_mesh(1,imaterial_id)
+ kappa_s = materials_ext_mesh(13,imaterial_id)
+ ! fluid properties
+ rho_f = materials_ext_mesh(2,imaterial_id)
+ kappa_f = materials_ext_mesh(14,imaterial_id)
+ eta_f = materials_ext_mesh(5,imaterial_id)
+ ! frame properties
+ kappa_fr = materials_ext_mesh(15,imaterial_id)
+ mu_fr = materials_ext_mesh(16,imaterial_id)
+ phi = materials_ext_mesh(3,imaterial_id)
+ tort = materials_ext_mesh(4,imaterial_id)
+ kxx = materials_ext_mesh(7,imaterial_id)
+ kxy = materials_ext_mesh(8,imaterial_id)
+ kxz = materials_ext_mesh(9,imaterial_id)
+ kyy = materials_ext_mesh(10,imaterial_id)
+ kyz = materials_ext_mesh(11,imaterial_id)
+ kzz = materials_ext_mesh(12,imaterial_id)
+ endif !if(idomain_id == IDOMAIN_ACOUSTIC .or. idomain_id == IDOMAIN_ELASTIC)
+
else if (mat_ext_mesh(2,ispec) == 1) then
stop 'material: interface not implemented yet'
@@ -196,6 +227,8 @@
! stores velocity model
+ if(idomain_id == IDOMAIN_ACOUSTIC .or. idomain_id == IDOMAIN_ELASTIC) then ! elastic or acoustic
+
! density
rhostore(i,j,k,ispec) = rho
@@ -209,8 +242,58 @@
! Stacey, a completer par la suite
rho_vp(i,j,k,ispec) = rho*vp
rho_vs(i,j,k,ispec) = rho*vs
+ !
+ rho_vpI(i,j,k,ispec) = rho*vp
+ rho_vpII(i,j,k,ispec) = 0.d0
+ rho_vsI(i,j,k,ispec) = rho*vs
+ rhoarraystore(1,i,j,k,ispec) = rho
+ rhoarraystore(2,i,j,k,ispec) = rho
+ phistore(i,j,k,ispec) = 0.d0
+ tortstore(i,j,k,ispec) = 1.d0
!end pll
+ else ! poroelastic
+
+ ! solid properties
+ rhoarraystore(1,i,j,k,ispec) = rho_s
+ kappaarraystore(1,i,j,k,ispec) = kappa_s
+ ! fluid properties
+ rhoarraystore(2,i,j,k,ispec) = rho_f
+ kappaarraystore(2,i,j,k,ispec) = kappa_f
+ etastore(i,j,k,ispec) = eta_f
+ ! frame properties
+ kappaarraystore(3,i,j,k,ispec) = kappa_fr
+ mustore(i,j,k,ispec) = mu_fr
+ phistore(i,j,k,ispec) = phi
+ tortstore(i,j,k,ispec) = tort
+ permstore(1,i,j,k,ispec) = kxx
+ permstore(2,i,j,k,ispec) = kxy
+ permstore(3,i,j,k,ispec) = kxz
+ permstore(4,i,j,k,ispec) = kyy
+ permstore(5,i,j,k,ispec) = kyz
+ permstore(6,i,j,k,ispec) = kzz
+
+ !Biot coefficients for the input phi
+ D_biot = kappa_s*(1._CUSTOM_REAL + phi*(kappa_s/kappa_f - 1._CUSTOM_REAL))
+ H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + 4._CUSTOM_REAL*mu_fr/3._CUSTOM_REAL
+ C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
+ M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
+ ! Approximated velocities (no viscous dissipation)
+ rho_bar = (1._CUSTOM_REAL - phi)*rho_s + phi*rho_f
+ afactor = rho_bar - phi/tort*rho_f
+ bfactor = H_biot + phi*rho_bar/(tort*rho_f)*M_biot - TWO*phi/tort*C_biot
+ cfactor = phi/(tort*rho_f)*(H_biot*M_biot - C_biot*C_biot)
+ cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(2._CUSTOM_REAL*afactor)
+ cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(2._CUSTOM_REAL*afactor)
+ cssquare = mu_fr/afactor
+
+ ! AC based on cpI,cpII & cs
+ rho_vpI(i,j,k,ispec) = (rho_bar - phi/tort*rho_f)*sqrt(cpIsquare)
+ rho_vpII(i,j,k,ispec) = (rho_bar - phi/tort*rho_f)*sqrt(cpIIsquare)
+ rho_vsI(i,j,k,ispec) = (rho_bar - phi/tort*rho_f)*sqrt(cssquare)
+
+ endif !if(idomain_id == IDOMAIN_ACOUSTIC .or. idomain_id == IDOMAIN_ELASTIC)
+
! adds anisotropic perturbation to vp, vs
if( ANISOTROPY ) then
!call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
@@ -252,7 +335,6 @@
else if( idomain_id == IDOMAIN_ELASTIC ) then
ispec_is_elastic(ispec) = .true.
else if( idomain_id == IDOMAIN_POROELASTIC ) then
- stop 'poroelastic material domain not implemented yet'
ispec_is_poroelastic(ispec) = .true.
else
stop 'error material domain index'
@@ -277,7 +359,11 @@
stop 'error material domain index element'
endif
! checks if domain is unique
- if( ispec_is_acoustic(ispec) .eqv. .true. .and. ispec_is_elastic(ispec) .eqv. .true. ) then
+ if( ((ispec_is_acoustic(ispec) .eqv. .true.) .and. (ispec_is_elastic(ispec) .eqv. .true.)) .or. &
+ ((ispec_is_acoustic(ispec) .eqv. .true.) .and. (ispec_is_poroelastic(ispec) .eqv. .true.)) .or. &
+ ((ispec_is_poroelastic(ispec) .eqv. .true.) .and. (ispec_is_elastic(ispec) .eqv. .true.)) .or. &
+ ((ispec_is_acoustic(ispec) .eqv. .true.) .and. (ispec_is_elastic(ispec) .eqv. .true.) .and. &
+ (ispec_is_poroelastic(ispec) .eqv. .true.)) ) then
print*,'error material domain assigned twice to element:',ispec
print*,'acoustic: ',ispec_is_acoustic(ispec)
print*,'elastic: ',ispec_is_elastic(ispec)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -32,6 +32,8 @@
gammaxstore,gammaystore,gammazstore, &
jacobianstore, rho_vp,rho_vs,qmu_attenuation_store, &
rhostore,kappastore,mustore, &
+ rhoarraystore,kappaarraystore,etastore,phistore,tortstore,permstore, &
+ rho_vpI,rho_vpII,rho_vsI, &
rmass,rmass_acoustic,rmass_solid_poroelastic,rmass_fluid_poroelastic, &
OCEANS,rmass_ocean_load,NGLOB_OCEAN,&
ibool, &
@@ -45,6 +47,10 @@
coupling_ac_el_normal,coupling_ac_el_jacobian2Dw, &
coupling_ac_el_ijk,coupling_ac_el_ispec, &
num_coupling_ac_el_faces, &
+ coupling_ac_po_normal,coupling_ac_po_jacobian2Dw, &
+ 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, &
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, &
@@ -71,6 +77,11 @@
! material
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rhostore,kappastore,mustore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: etastore,phistore,tortstore
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,nspec) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,nspec) :: kappaarraystore
+ real(kind=CUSTOM_REAL), dimension(6,NGLLX,NGLLY,NGLLZ,nspec) :: permstore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rho_vpI,rho_vpII,rho_vsI
real(kind=CUSTOM_REAL), dimension(nglob) :: rmass,rmass_acoustic, &
rmass_solid_poroelastic,rmass_fluid_poroelastic
! ocean load
@@ -103,6 +114,20 @@
integer :: coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces)
integer :: coupling_ac_el_ispec(num_coupling_ac_el_faces)
+! acoustic-poroelastic coupling surface
+ integer :: num_coupling_ac_po_faces
+ real(kind=CUSTOM_REAL) :: coupling_ac_po_normal(NDIM,NGLLSQUARE,num_coupling_ac_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_ac_po_jacobian2Dw(NGLLSQUARE,num_coupling_ac_po_faces)
+ integer :: coupling_ac_po_ijk(3,NGLLSQUARE,num_coupling_ac_po_faces)
+ integer :: coupling_ac_po_ispec(num_coupling_ac_po_faces)
+
+! 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_el_po_ispec(num_coupling_el_po_faces)
+
! MPI interfaces
integer :: num_interfaces_ext_mesh
integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
@@ -197,6 +222,15 @@
if( POROELASTIC_SIMULATION ) then
write(IOUT) rmass_solid_poroelastic
write(IOUT) rmass_fluid_poroelastic
+ write(IOUT) rhoarraystore
+ write(IOUT) kappaarraystore
+ write(IOUT) etastore
+ write(IOUT) tortstore
+ write(IOUT) permstore
+ write(IOUT) phistore
+ write(IOUT) rho_vpI
+ write(IOUT) rho_vpII
+ write(IOUT) rho_vsI
endif
! absorbing boundary surface
@@ -220,6 +254,20 @@
write(IOUT) coupling_ac_el_jacobian2Dw
write(IOUT) coupling_ac_el_normal
+! acoustic-poroelastic coupling surface
+ write(IOUT) num_coupling_ac_po_faces
+ write(IOUT) coupling_ac_po_ispec
+ write(IOUT) coupling_ac_po_ijk
+ write(IOUT) coupling_ac_po_jacobian2Dw
+ write(IOUT) coupling_ac_po_normal
+
+! elastic-poroelastic coupling surface
+ write(IOUT) num_coupling_el_po_faces
+ write(IOUT) coupling_el_po_ispec
+ write(IOUT) coupling_el_po_ijk
+ write(IOUT) coupling_el_po_jacobian2Dw
+ write(IOUT) coupling_el_po_normal
+
!MPI interfaces
write(IOUT) num_interfaces_ext_mesh
write(IOUT) maxval(nibool_interfaces_ext_mesh(:))
@@ -418,7 +466,90 @@
deallocate(iglob_tmp)
endif
+ ! acoustic-poroelastic domains
+ if( ACOUSTIC_SIMULATION .and. POROELASTIC_SIMULATION ) then
+ ! saves points on acoustic-poroelastic coupling interface
+ allocate( iglob_tmp(NGLLSQUARE*num_coupling_ac_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_tmp'
+ inum = 0
+ iglob_tmp(:) = 0
+ do i=1,num_coupling_ac_po_faces
+ do j=1,NGLLSQUARE
+ inum = inum+1
+ iglob_tmp(inum) = ibool(coupling_ac_po_ijk(1,j,i), &
+ coupling_ac_po_ijk(2,j,i), &
+ coupling_ac_po_ijk(3,j,i), &
+ coupling_ac_po_ispec(i) )
+ enddo
+ enddo
+ filename = prname(1:len_trim(prname))//'coupling_acoustic_poroelastic'
+ call write_VTK_data_points(nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iglob_tmp,NGLLSQUARE*num_coupling_ac_po_faces, &
+ filename)
+ ! saves acoustic/poroelastic flag
+ allocate(v_tmp_i(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array v_tmp_i'
+ do i=1,nspec
+ if( ispec_is_acoustic(i) ) then
+ v_tmp_i(i) = 1
+ else if( ispec_is_poroelastic(i) ) then
+ v_tmp_i(i) = 2
+ else
+ v_tmp_i(i) = 0
+ endif
+ enddo
+ filename = prname(1:len_trim(prname))//'acoustic_poroelastic_flag'
+ call write_VTK_data_elem_i(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ v_tmp_i,filename)
+
+ deallocate(v_tmp_i,iglob_tmp)
+ endif !if( ACOUSTIC_SIMULATION .and. POROELASTIC_SIMULATION )
+
+ ! elastic-poroelastic domains
+ if( ELASTIC_SIMULATION .and. POROELASTIC_SIMULATION ) then
+ ! saves points on elastic-poroelastic coupling interface
+ allocate( iglob_tmp(NGLLSQUARE*num_coupling_el_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_tmp'
+ inum = 0
+ iglob_tmp(:) = 0
+ do i=1,num_coupling_el_po_faces
+ do j=1,NGLLSQUARE
+ inum = inum+1
+ iglob_tmp(inum) = ibool(coupling_el_po_ijk(1,j,i), &
+ coupling_el_po_ijk(2,j,i), &
+ coupling_el_po_ijk(3,j,i), &
+ coupling_el_po_ispec(i) )
+ enddo
+ enddo
+ filename = prname(1:len_trim(prname))//'coupling_elastic_poroelastic'
+ call write_VTK_data_points(nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iglob_tmp,NGLLSQUARE*num_coupling_el_po_faces, &
+ filename)
+
+ ! saves elastic/poroelastic flag
+ allocate(v_tmp_i(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array v_tmp_i'
+ do i=1,nspec
+ if( ispec_is_elastic(i) ) then
+ v_tmp_i(i) = 1
+ else if( ispec_is_poroelastic(i) ) then
+ v_tmp_i(i) = 2
+ else
+ v_tmp_i(i) = 0
+ endif
+ enddo
+ filename = prname(1:len_trim(prname))//'elastic_poroelastic_flag'
+ call write_VTK_data_elem_i(nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
+ v_tmp_i,filename)
+
+ deallocate(v_tmp_i,iglob_tmp)
+ endif !if( ACOUSTIC_SIMULATION .and. POROELASTIC_SIMULATION
+
!! saves 1. MPI interface
! if( num_interfaces_ext_mesh >= 1 ) then
! filename = prname(1:len_trim(prname))//'MPI_1_points'
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -27,7 +27,8 @@
subroutine check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
kappastore,mustore,rho_vp,rho_vs, &
- DT, model_speed_max,min_resolved_period )
+ DT, model_speed_max,min_resolved_period, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI )
! check the mesh, stability and resolved period
!
@@ -39,6 +40,9 @@
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) :: 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
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
double precision :: DT
@@ -46,6 +50,7 @@
! local parameters
real(kind=CUSTOM_REAL) :: vpmin,vpmax,vsmin,vsmax,vpmin_glob,vpmax_glob,vsmin_glob,vsmax_glob
+ real(kind=CUSTOM_REAL) :: vp2min,vp2max,vp2min_glob,vp2max_glob
real(kind=CUSTOM_REAL) :: distance_min,distance_max,distance_min_glob,distance_max_glob
real(kind=CUSTOM_REAL) :: elemsize_min,elemsize_max,elemsize_min_glob,elemsize_max_glob
real(kind=CUSTOM_REAL) :: cmax,cmax_glob,pmax,pmax_glob
@@ -70,7 +75,7 @@
!real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
- logical :: has_vs_zero
+ logical :: has_vs_zero,has_vp2_zero
! initializations
if( DT <= 0.0d0) then
@@ -82,6 +87,9 @@
vpmin_glob = HUGEVAL
vpmax_glob = -HUGEVAL
+ vp2min_glob = HUGEVAL
+ vp2max_glob = -HUGEVAL
+
vsmin_glob = HUGEVAL
vsmax_glob = -HUGEVAL
@@ -97,18 +105,23 @@
dt_suggested_glob = HUGEVAL
has_vs_zero = .false.
+ has_vp2_zero = .false.
! checks courant number & minimum resolved period for each grid cell
do ispec=1,NSPEC_AB
! determines minimum/maximum velocities within this element
- call get_vpvs_minmax(vpmin,vpmax,vsmin,vsmax,ispec,has_vs_zero, &
- NSPEC_AB,kappastore,mustore,rho_vp,rho_vs)
+ call get_vpvs_minmax(vpmin,vpmax,vp2min,vp2max,vsmin,vsmax,ispec,has_vs_zero, &
+ has_vp2_zero,NSPEC_AB,kappastore,mustore,rho_vp,rho_vs, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI)
! min/max for whole cpu partition
vpmin_glob = min ( vpmin_glob, vpmin)
vpmax_glob = max ( vpmax_glob, vpmax)
+ vp2min_glob = min ( vp2min_glob, vp2min)
+ vp2max_glob = max ( vp2max_glob, vp2max)
+
vsmin_glob = min ( vsmin_glob, vsmin)
vsmax_glob = max ( vsmax_glob, vsmax)
@@ -130,12 +143,12 @@
! based on minimum GLL point distance and maximum velocity
! i.e. on the maximum ratio of ( velocity / gridsize )
if( DT_PRESENT ) then
- cmax = max( vpmax,vsmax ) * DT / distance_min
+ cmax = max( vpmax,vp2max,vsmax ) * DT / distance_min
cmax_glob = max(cmax_glob,cmax)
endif
! suggested timestep
- dt_suggested = COURANT_SUGGESTED * distance_min / max( vpmax,vsmax )
+ dt_suggested = COURANT_SUGGESTED * distance_min / max( vpmax,vp2max,vsmax )
dt_suggested_glob = min( dt_suggested_glob, dt_suggested)
! estimation of minimum period resolved
@@ -148,14 +161,14 @@
avg_distance = elemsize_max / ( NGLLX - 1 ) ! since NGLLX = NGLLY = NGLLZ
! biggest possible minimum period such that number of points per minimum wavelength
- ! npts = ( min(vpmin,vsmin) * pmax ) / avg_distance is about ~ NPTS_PER_WAVELENGTH
+ ! npts = ( min(vpmin,vp2min,vsmin) * pmax ) / avg_distance is about ~ NPTS_PER_WAVELENGTH
!
! note: obviously, this estimation depends on the choice of points per wavelength
! which is empirical at the moment.
! also, keep in mind that the minimum period is just an estimation and
! there is no such sharp cut-off period for valid synthetics.
! seismograms become just more and more inaccurate for periods shorter than this estimate.
- pmax = avg_distance / min( vpmin,vsmin ) * NPTS_PER_WAVELENGTH
+ pmax = avg_distance / min( vpmin,vp2max,vsmin ) * NPTS_PER_WAVELENGTH
pmax_glob = max(pmax_glob,pmax)
@@ -186,6 +199,14 @@
call min_all_cr(vpmin,vpmin_glob)
call max_all_cr(vpmax,vpmax_glob)
+ ! Vp2 velocity (relevant for poroelastic cases)
+ vp2min = vp2min_glob
+ if( has_vp2_zero ) vp2min = 0.0
+
+ vp2max = vp2max_glob
+ call min_all_cr(vp2min,vp2min_glob)
+ call max_all_cr(vp2max,vp2max_glob)
+
! Vs velocity
vsmin = vsmin_glob
if( has_vs_zero ) vsmin = 0.0
@@ -201,6 +222,12 @@
if( vpmax_glob >= HUGEVAL ) then
call exit_mpi(myrank,"error: vp maximum velocity")
endif
+ if( vp2min_glob < 0.0_CUSTOM_REAL ) then
+ call exit_mpi(myrank,"error: vp2 minimum velocity")
+ endif
+ if( vp2max_glob >= HUGEVAL ) then
+ call exit_mpi(myrank,"error: vp2 maximum velocity")
+ endif
if( vsmin_glob < 0.0_CUSTOM_REAL ) then
call exit_mpi(myrank,"error: vs minimum velocity")
endif
@@ -264,6 +291,7 @@
write(IMAIN,*)
write(IMAIN,*) '********'
write(IMAIN,*) 'Model: P velocity min,max = ',vpmin_glob,vpmax_glob
+ write(IMAIN,*) 'Model: PII velocity min,max = ',vp2min_glob,vp2max_glob
write(IMAIN,*) 'Model: S velocity min,max = ',vsmin_glob,vsmax_glob
write(IMAIN,*) '********'
write(IMAIN,*)
@@ -310,31 +338,37 @@
!
- subroutine get_vpvs_minmax(vpmin,vpmax,vsmin,vsmax,ispec,has_vs_zero, &
- NSPEC_AB,kappastore,mustore,rho_vp,rho_vs)
-
+ subroutine get_vpvs_minmax(vpmin,vpmax,vp2min,vp2max,vsmin,vsmax,ispec,has_vs_zero, &
+ has_vp2_zero,NSPEC_AB,kappastore,mustore,rho_vp,rho_vs, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI)
! calculates the min/max size of the specified element (ispec)
implicit none
include "constants.h"
- real(kind=CUSTOM_REAL) :: vpmin,vpmax,vsmin,vsmax
+ real(kind=CUSTOM_REAL) :: vpmin,vpmax,vp2min,vp2max,vsmin,vsmax
integer :: ispec
logical :: has_vs_zero
+ logical :: has_vp2_zero
integer :: NSPEC_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) :: 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
! local parameters
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ):: vp_elem,vs_elem
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ):: vp_elem,vs_elem,vp2_elem
real(kind=CUSTOM_REAL), dimension(1) :: val_min,val_max
! initializes
vpmin = HUGEVAL
vpmax = -HUGEVAL
+ vp2min = HUGEVAL
+ vp2max = -HUGEVAL
vsmin = HUGEVAL
vsmax = -HUGEVAL
@@ -342,6 +376,10 @@
where( rho_vp(:,:,:,ispec) > TINYVAL )
vp_elem(:,:,:) = (FOUR_THIRDS * mustore(:,:,:,ispec) &
+ kappastore(:,:,:,ispec)) / rho_vp(:,:,:,ispec)
+ elsewhere(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) )
elsewhere
vp_elem(:,:,:) = 0.0
endwhere
@@ -352,9 +390,33 @@
vpmin = min(vpmin,val_min(1))
vpmax = max(vpmax,val_max(1))
+ ! vpII
+ where( rho_vpII(:,:,:,ispec) > TINYVAL)
+ vp2_elem(:,:,:) = rho_vpII(:,:,:,ispec)/( ((1._CUSTOM_REAL - phistore(:,:,:,ispec))* &
+ rhoarraystore(1,:,:,:,ispec) + phistore(:,:,:,ispec)*rhoarraystore(2,:,:,:,ispec)) - &
+ phistore(:,:,:,ispec)/tortstore(:,:,:,ispec)* rhoarraystore(2,:,:,:,ispec) )
+ elsewhere
+ vp2_elem(:,:,:) = 0.0
+ endwhere
+
+ val_min = minval(vp2_elem(:,:,:))
+ val_max = maxval(vp2_elem(:,:,:))
+
+ ! ignore non porous region with vp2 = 0
+ if( val_min(1) > 0.0001 ) then
+ vp2min = min(vp2min,val_min(1))
+ else
+ has_vp2_zero = .true.
+ endif
+ vp2max = max(vp2max,val_max(1))
+
! vs
where( rho_vs(:,:,:,ispec) > TINYVAL )
vs_elem(:,:,:) = mustore(:,:,:,ispec) / rho_vs(:,:,:,ispec)
+ elsewhere(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) )
elsewhere
vs_elem(:,:,:) = 0.0
endwhere
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -255,3 +255,152 @@
endif
end subroutine assemble_MPI_vector_ext_mesh_w
+
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine assemble_MPI_vector_ext_mesh_po_s(NPROC,NGLOB_AB,array_val1,array_val2, &
+ buffer_send_vector_ext_mesh_s,buffer_recv_vector_ext_mesh_s, &
+ buffer_send_vector_ext_mesh_w,buffer_recv_vector_ext_mesh_w, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh_s,request_recv_vector_ext_mesh_s, &
+ request_send_vector_ext_mesh_w,request_recv_vector_ext_mesh_w &
+ )
+
+! sends data
+
+ implicit none
+
+ include "constants.h"
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val1,array_val2
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
+ buffer_send_vector_ext_mesh_s,buffer_recv_vector_ext_mesh_s,&
+ buffer_send_vector_ext_mesh_w,buffer_recv_vector_ext_mesh_w
+
+ integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh_s,request_recv_vector_ext_mesh_s
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh_w,request_recv_vector_ext_mesh_w
+
+ integer ipoin,iinterface
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! partition border copy into the buffer
+ do iinterface = 1, num_interfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ buffer_send_vector_ext_mesh_s(:,ipoin,iinterface) = array_val1(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
+ buffer_send_vector_ext_mesh_w(:,ipoin,iinterface) = array_val2(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
+ enddo
+ enddo
+
+! send messages
+ do iinterface = 1, num_interfaces_ext_mesh
+!val1
+ call isend_cr(buffer_send_vector_ext_mesh_s(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_vector_ext_mesh_s(iinterface) &
+ )
+ call irecv_cr(buffer_recv_vector_ext_mesh_s(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_vector_ext_mesh_s(iinterface) &
+ )
+!val2
+ call isend_cr(buffer_send_vector_ext_mesh_w(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_vector_ext_mesh_w(iinterface) &
+ )
+ call irecv_cr(buffer_recv_vector_ext_mesh_w(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_vector_ext_mesh_w(iinterface) &
+ )
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_ext_mesh_po_s
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine assemble_MPI_vector_ext_mesh_po_w(NPROC,NGLOB_AB,array_val1,array_val2, &
+ buffer_recv_vector_ext_mesh_s,buffer_recv_vector_ext_mesh_w, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ request_send_vector_ext_mesh_s,request_recv_vector_ext_mesh_s, &
+ request_send_vector_ext_mesh_w,request_recv_vector_ext_mesh_w)
+
+! waits for data to receive and assembles
+
+ implicit none
+
+ include "constants.h"
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val1,array_val2
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: &
+ buffer_recv_vector_ext_mesh_s,buffer_recv_vector_ext_mesh_w
+
+ integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh_s,request_recv_vector_ext_mesh_s
+ integer, dimension(num_interfaces_ext_mesh) :: request_send_vector_ext_mesh_w,request_recv_vector_ext_mesh_w
+
+ integer ipoin,iinterface
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! wait for communications completion (recv)
+ do iinterface = 1, num_interfaces_ext_mesh
+ call wait_req(request_recv_vector_ext_mesh_s(iinterface))
+ call wait_req(request_recv_vector_ext_mesh_w(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ do iinterface = 1, num_interfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ array_val1(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val1(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh_s(:,ipoin,iinterface)
+ array_val2(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val2(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh_w(:,ipoin,iinterface)
+ enddo
+ enddo
+
+! wait for communications completion (send)
+ do iinterface = 1, num_interfaces_ext_mesh
+ call wait_req(request_send_vector_ext_mesh_s(iinterface))
+ call wait_req(request_send_vector_ext_mesh_w(iinterface))
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_ext_mesh_po_w
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -81,7 +81,7 @@
if( ELASTIC_SIMULATION ) call compute_forces_elastic()
! poroelastic solver
- if( POROELASTIC_SIMULATION ) stop 'poroelastic simulation not implemented yet'
+ if( POROELASTIC_SIMULATION ) call compute_forces_poroelastic()
! restores last time snapshot saved for backward/reconstruction of wavefields
! note: this must be read in after the Newark time scheme
@@ -136,6 +136,7 @@
use specfem_par
use specfem_par_elastic
+ use specfem_par_poroelastic
use specfem_par_acoustic
implicit none
@@ -144,17 +145,27 @@
ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
ihours_total,iminutes_total,iseconds_total,int_t_total
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!chris: Rewrite to get norm for each material when coupled simulations
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
! compute maximum of norm of displacement in each slice
- if( ELASTIC_SIMULATION ) then
+ if( ELASTIC_SIMULATION ) &
Usolidnorm = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2 + displ(3,:)**2))
- else
- if( ACOUSTIC_SIMULATION ) then
- Usolidnorm = maxval(abs(potential_dot_dot_acoustic(:)))
- endif
+ if( ACOUSTIC_SIMULATION ) &
+ Usolidnormp = maxval(abs(potential_dot_dot_acoustic(:)))
+ if( POROELASTIC_SIMULATION ) then
+ Usolidnorms = maxval(sqrt(displs_poroelastic(1,:)**2 + displs_poroelastic(2,:)**2 + &
+ displs_poroelastic(3,:)**2))
+ Usolidnormw = maxval(sqrt(displw_poroelastic(1,:)**2 + displw_poroelastic(2,:)**2 + &
+ displw_poroelastic(3,:)**2))
endif
! compute the maximum of the maxima for all the slices using an MPI reduction
call max_all_cr(Usolidnorm,Usolidnorm_all)
+ call max_all_cr(Usolidnormp,Usolidnormp_all)
+ call max_all_cr(Usolidnorms,Usolidnorms_all)
+ call max_all_cr(Usolidnormw,Usolidnormw_all)
! adjoint simulations
if( SIMULATION_TYPE == 3 ) then
@@ -183,12 +194,13 @@
write(IMAIN,*) 'Elapsed time in seconds = ',tCPU
write(IMAIN,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
- if( ELASTIC_SIMULATION ) then
+ if( ELASTIC_SIMULATION ) &
write(IMAIN,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
- else
- if( ACOUSTIC_SIMULATION ) then
- write(IMAIN,*) 'Max norm pressure P in all slices (Pa) = ',Usolidnorm_all
- endif
+ if( ACOUSTIC_SIMULATION ) &
+ write(IMAIN,*) 'Max norm pressure P in all slices (Pa) = ',Usolidnormp_all
+ if( POROELASTIC_SIMULATION ) then
+ write(IMAIN,*) 'Max norm displacement vector Us in all slices (m) = ',Usolidnorms_all
+ write(IMAIN,*) 'Max norm displacement vector W in all slices (m) = ',Usolidnormw_all
endif
! adjoint simulations
if (SIMULATION_TYPE == 3) write(IMAIN,*) &
@@ -233,7 +245,14 @@
write(IOUT,*) 'Elapsed time in seconds = ',tCPU
write(IOUT,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
- write(IOUT,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
+ if( ELASTIC_SIMULATION ) &
+ write(IOUT,*) 'Max norm displacement vector U in all slices (m) = ',Usolidnorm_all
+ if( ACOUSTIC_SIMULATION ) &
+ write(IOUT,*) 'Max norm pressure P in all slices (Pa) = ',Usolidnormp_all
+ if( POROELASTIC_SIMULATION ) then
+ write(IOUT,*) 'Max norm displacement vector Us in all slices (m) = ',Usolidnorms_all
+ write(IOUT,*) 'Max norm displacement vector W in all slices (m) = ',Usolidnormw_all
+ endif
! adjoint simulations
if (SIMULATION_TYPE == 3) write(IOUT,*) &
'Max norm U (backward) in all slices = ',b_Usolidnorm_all
@@ -243,7 +262,10 @@
! check stability of the code, exit if unstable
! negative values can occur with some compilers when the unstable value is greater
! than the greatest possible floating-point number of the machine
- if(Usolidnorm_all > STABILITY_THRESHOLD .or. Usolidnorm_all < 0) &
+ if(Usolidnorm_all > STABILITY_THRESHOLD .or. Usolidnorm_all < 0 &
+ .or. Usolidnormp_all > STABILITY_THRESHOLD .or. Usolidnormp_all < 0 &
+ .or. Usolidnorms_all > STABILITY_THRESHOLD .or. Usolidnorms_all < 0 &
+ .or. Usolidnormw_all > STABILITY_THRESHOLD .or. Usolidnormw_all < 0) &
call exit_MPI(myrank,'forward simulation became unstable and blew up')
! adjoint simulations
if(SIMULATION_TYPE == 3 .and. (b_Usolidnorm_all > STABILITY_THRESHOLD &
@@ -328,6 +350,21 @@
accel(:,:) = 0._CUSTOM_REAL
endif
+! updates poroelastic displacements and velocities
+ if( POROELASTIC_SIMULATION ) then
+ ! solid phase
+ displs_poroelastic(:,:) = displs_poroelastic(:,:) + deltat*velocs_poroelastic(:,:) + &
+ deltatsqover2*accels_poroelastic(:,:)
+ velocs_poroelastic(:,:) = velocs_poroelastic(:,:) + deltatover2*accels_poroelastic(:,:)
+ accels_poroelastic(:,:) = 0._CUSTOM_REAL
+
+ ! fluid phase
+ displw_poroelastic(:,:) = displw_poroelastic(:,:) + deltat*velocw_poroelastic(:,:) + &
+ deltatsqover2*accelw_poroelastic(:,:)
+ velocw_poroelastic(:,:) = velocw_poroelastic(:,:) + deltatover2*accelw_poroelastic(:,:)
+ accelw_poroelastic(:,:) = 0._CUSTOM_REAL
+ endif
+
! adjoint simulations
if (SIMULATION_TYPE == 3) then
! acoustic backward fields
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -37,7 +37,7 @@
UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
PRINT_SOURCE_TIME_FUNCTION, &
nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh, &
- ispec_is_acoustic,ispec_is_elastic, &
+ ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
implicit none
@@ -58,7 +58,7 @@
! arrays containing coordinates of the points
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
- logical, dimension(NSPEC_AB) :: ispec_is_acoustic,ispec_is_elastic
+ logical, dimension(NSPEC_AB) :: ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic
integer yr,jda,ho,mi
@@ -405,6 +405,8 @@
idomain(isource) = IDOMAIN_ACOUSTIC
else if( ispec_is_elastic( ispec_selected_source(isource) ) ) then
idomain(isource) = IDOMAIN_ELASTIC
+ else if( ispec_is_poroelastic( ispec_selected_source(isource) ) ) then
+ idomain(isource) = IDOMAIN_POROELASTIC
else
idomain(isource) = 0
endif
@@ -776,6 +778,8 @@
write(IMAIN,*) ' in acoustic domain'
else if( idomain(isource) == IDOMAIN_ELASTIC ) then
write(IMAIN,*) ' in elastic domain'
+ else if( idomain(isource) == IDOMAIN_POROELASTIC ) then
+ write(IMAIN,*) ' in poroelastic domain'
else
write(IMAIN,*) ' in unknown domain'
endif
@@ -887,8 +891,8 @@
endif
! checks source domain
- if( idomain(isource) /= IDOMAIN_ACOUSTIC .and. idomain(isource) /= IDOMAIN_ELASTIC ) then
- ! only acoustic/elastic domain implement yet
+ if( idomain(isource) /= IDOMAIN_ACOUSTIC .and. idomain(isource) /= IDOMAIN_ELASTIC .and. &
+ idomain(isource) /= IDOMAIN_POROELASTIC ) then
call exit_MPI(myrank,'source located in unknown domain')
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -132,6 +132,18 @@
if(FIX_UNDERFLOW_PROBLEM) displ(:,:) = VERYSMALLVAL
endif
+ ! initialize poroelastic arrays to zero/verysmallvall
+ if( POROELASTIC_SIMULATION ) then
+ displs_poroelastic(:,:) = 0._CUSTOM_REAL
+ velocs_poroelastic(:,:) = 0._CUSTOM_REAL
+ accels_poroelastic(:,:) = 0._CUSTOM_REAL
+ displw_poroelastic(:,:) = 0._CUSTOM_REAL
+ velocw_poroelastic(:,:) = 0._CUSTOM_REAL
+ accelw_poroelastic(:,:) = 0._CUSTOM_REAL
+ ! put negligible initial value to avoid very slow underflow trapping
+ if(FIX_UNDERFLOW_PROBLEM) displs_poroelastic(:,:) = VERYSMALLVAL
+ if(FIX_UNDERFLOW_PROBLEM) displw_poroelastic(:,:) = VERYSMALLVAL
+ endif
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
@@ -251,9 +263,6 @@
endif ! ELASTIC_SIMULATION
if(POROELASTIC_SIMULATION) then
-
- stop 'poroelastic simulation not implemented yet'
- ! but would be something like this...
call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_solid_poroelastic, &
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -187,16 +187,47 @@
! poroelastic
call any_all_l( ANY(ispec_is_poroelastic), POROELASTIC_SIMULATION )
if( POROELASTIC_SIMULATION ) then
+ ! displacement,velocity,acceleration for the solid (s) & fluid (w) phases
+ allocate(displs_poroelastic(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array displs_poroelastic'
+ allocate(velocs_poroelastic(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array velocs_poroelastic'
+ allocate(accels_poroelastic(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array accels_poroelastic'
+ allocate(displw_poroelastic(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array displw_poroelastic'
+ allocate(velocw_poroelastic(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array velocw_poroelastic'
+ allocate(accelw_poroelastic(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array accelw_poroelastic'
- stop 'not implemented yet: read rmass_solid_poroelastic .. '
-
allocate(rmass_solid_poroelastic(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array rmass_solid_poroelastic'
allocate(rmass_fluid_poroelastic(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array rmass_fluid_poroelastic'
+ allocate(rhoarraystore(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ kappaarraystore(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ etastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ tortstore(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ phistore(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ permstore(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ rho_vpI(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ rho_vpII(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ rho_vsI(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array poroelastic properties'
+
read(27) rmass_solid_poroelastic
read(27) rmass_fluid_poroelastic
+ read(27) rhoarraystore
+ read(27) kappaarraystore
+ read(27) etastore
+ read(27) tortstore
+ read(27) permstore
+ read(27) phistore
+ read(27) rho_vpI
+ read(27) rho_vpII
+ read(27) rho_vsI
endif
! checks simulation types are valid
@@ -243,6 +274,30 @@
read(27) coupling_ac_el_jacobian2Dw
read(27) coupling_ac_el_normal
+! acoustic-poroelastic coupling surface
+ read(27) num_coupling_ac_po_faces
+ allocate(coupling_ac_po_normal(NDIM,NGLLSQUARE,num_coupling_ac_po_faces), &
+ coupling_ac_po_jacobian2Dw(NGLLSQUARE,num_coupling_ac_po_faces), &
+ coupling_ac_po_ijk(3,NGLLSQUARE,num_coupling_ac_po_faces), &
+ coupling_ac_po_ispec(num_coupling_ac_po_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_po_normal etc.'
+ read(27) coupling_ac_po_ispec
+ read(27) coupling_ac_po_ijk
+ read(27) coupling_ac_po_jacobian2Dw
+ read(27) coupling_ac_po_normal
+
+! elastic-poroelastic coupling surface
+ read(27) num_coupling_el_po_faces
+ 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)
+ if( ier /= 0 ) stop 'error allocating array coupling_el_po_normal etc.'
+ read(27) coupling_el_po_ispec
+ read(27) coupling_el_po_ijk
+ read(27) coupling_el_po_jacobian2Dw
+ read(27) coupling_el_po_normal
+
! MPI interfaces
read(27) num_interfaces_ext_mesh
read(27) max_nibool_interfaces_ext_mesh
@@ -288,7 +343,15 @@
request_send_vector_ext_mesh(num_interfaces_ext_mesh), &
request_recv_vector_ext_mesh(num_interfaces_ext_mesh), &
request_send_scalar_ext_mesh(num_interfaces_ext_mesh), &
- request_recv_scalar_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ request_recv_scalar_ext_mesh(num_interfaces_ext_mesh), &
+ buffer_send_vector_ext_mesh_s(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ buffer_recv_vector_ext_mesh_s(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ buffer_send_vector_ext_mesh_w(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ buffer_recv_vector_ext_mesh_w(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ request_send_vector_ext_mesh_s(num_interfaces_ext_mesh), &
+ request_recv_vector_ext_mesh_s(num_interfaces_ext_mesh), &
+ request_send_vector_ext_mesh_w(num_interfaces_ext_mesh), &
+ request_recv_vector_ext_mesh_w(num_interfaces_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating array buffer_send_vector_ext_mesh etc.'
! locate inner and outer elements
@@ -311,8 +374,25 @@
! check courant criteria on mesh
if( ELASTIC_SIMULATION ) then
+ allocate(rho_vpI(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rho_vpII(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rho_vsI(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ rho_vpI = 0.0_CUSTOM_REAL
+ rho_vpII = 0.0_CUSTOM_REAL
+ rho_vsI = 0.0_CUSTOM_REAL
call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
- kappastore,mustore,rho_vp,rho_vs,DT,model_speed_max,min_resolved_period )
+ kappastore,mustore,rho_vp,rho_vs,DT,model_speed_max,min_resolved_period, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI )
+ deallocate(rho_vpI,rho_vpII,rho_vsI)
+ else if( POROELASTIC_SIMULATION ) then
+ allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ rho_vp = 0.0_CUSTOM_REAL
+ rho_vs = 0.0_CUSTOM_REAL
+ call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
+ kappastore,mustore,rho_vp,rho_vs,DT,model_speed_max,min_resolved_period, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI )
+ deallocate(rho_vp,rho_vs)
else if( ACOUSTIC_SIMULATION ) then
allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array rho_vp'
@@ -321,7 +401,8 @@
rho_vp = sqrt( kappastore / rhostore ) * rhostore
rho_vs = 0.0_CUSTOM_REAL
call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
- kappastore,mustore,rho_vp,rho_vs,DT,model_speed_max,min_resolved_period )
+ kappastore,mustore,rho_vp,rho_vs,DT,model_speed_max,min_resolved_period, &
+ phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI )
deallocate(rho_vp,rho_vs)
endif
@@ -337,6 +418,7 @@
use specfem_par
use specfem_par_elastic
+ use specfem_par_poroelastic
use specfem_par_acoustic
implicit none
! local parameters
@@ -455,6 +537,42 @@
!print *,'rank ',myrank,' elastic outer spec: ',nspec_outer_elastic
endif
+! sets up elements for loops in poroelastic simulations
+ if( POROELASTIC_SIMULATION ) then
+ ! counts inner and outer elements
+ nspec_inner_poroelastic = 0
+ nspec_outer_poroelastic = 0
+ do ispec = 1, NSPEC_AB
+ if( ispec_is_poroelastic(ispec) ) then
+ if( ispec_is_inner(ispec) .eqv. .true. ) then
+ nspec_inner_poroelastic = nspec_inner_poroelastic + 1
+ else
+ nspec_outer_poroelastic = nspec_outer_poroelastic + 1
+ endif
+ endif
+ enddo
+
+ ! stores indices of inner and outer elements for faster(?) computation
+ num_phase_ispec_poroelastic = max(nspec_inner_poroelastic,nspec_outer_poroelastic)
+ allocate( phase_ispec_inner_poroelastic(num_phase_ispec_poroelastic,2),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array phase_ispec_inner_poroelastic'
+ nspec_inner_poroelastic = 0
+ nspec_outer_poroelastic = 0
+ do ispec = 1, NSPEC_AB
+ if( ispec_is_poroelastic(ispec) ) then
+ if( ispec_is_inner(ispec) .eqv. .true. ) then
+ nspec_inner_poroelastic = nspec_inner_poroelastic + 1
+ phase_ispec_inner_poroelastic(nspec_inner_poroelastic,2) = ispec
+ else
+ nspec_outer_poroelastic = nspec_outer_poroelastic + 1
+ phase_ispec_inner_poroelastic(nspec_outer_poroelastic,1) = ispec
+ endif
+ endif
+ enddo
+ !print *,'rank ',myrank,' poroelastic inner spec: ',nspec_inner_poroelastic
+ !print *,'rank ',myrank,' poroelastic outer spec: ',nspec_outer_poroelastic
+ endif
+
end subroutine rmd_setup_inner_outer_elemnts
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -67,6 +67,7 @@
use specfem_par
use specfem_par_acoustic
use specfem_par_elastic
+ use specfem_par_poroelastic
use specfem_par_movie
implicit none
@@ -107,7 +108,7 @@
UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
PRINT_SOURCE_TIME_FUNCTION, &
nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh,&
- ispec_is_acoustic,ispec_is_elastic, &
+ ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
if(abs(minval(tshift_cmt)) > TINYVAL) call exit_MPI(myrank,'one tshift_cmt must be zero, others must be positive')
@@ -507,6 +508,7 @@
use specfem_par
use specfem_par_elastic
+ use specfem_par_poroelastic
use specfem_par_acoustic
implicit none
@@ -536,8 +538,8 @@
ispec = ispec_selected_source(isource)
- ! elastic moment tensor source
- if( ispec_is_elastic(ispec) ) then
+ ! elastic or poroelastic moment tensor source
+ if( ispec_is_elastic(ispec) .or. ispec_is_poroelastic(ispec)) then
call compute_arrays_source(ispec, &
xi_source(isource),eta_source(isource),gamma_source(isource),sourcearray, &
Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource), &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -179,6 +179,14 @@
integer, dimension(:), allocatable :: request_recv_scalar_ext_mesh
integer, dimension(:), allocatable :: request_send_vector_ext_mesh
integer, dimension(:), allocatable :: request_recv_vector_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_send_vector_ext_mesh_s
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_recv_vector_ext_mesh_s
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_send_vector_ext_mesh_w
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_recv_vector_ext_mesh_w
+ integer, dimension(:), allocatable :: request_send_vector_ext_mesh_s
+ integer, dimension(:), allocatable :: request_recv_vector_ext_mesh_s
+ integer, dimension(:), allocatable :: request_send_vector_ext_mesh_w
+ integer, dimension(:), allocatable :: request_recv_vector_ext_mesh_w
! for detecting surface receivers and source in case of external mesh
logical, dimension(:), allocatable :: iglob_is_surface_external_mesh
@@ -188,7 +196,10 @@
logical, dimension(:), allocatable :: ispec_is_inner
! maximum of the norm of the displacement
- real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all
+ real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all ! elastic
+ real(kind=CUSTOM_REAL) Usolidnormp,Usolidnormp_all ! acoustic
+ real(kind=CUSTOM_REAL) Usolidnorms,Usolidnorms_all ! solid poroelastic
+ real(kind=CUSTOM_REAL) Usolidnormw,Usolidnormw_all ! fluid (w.r.t.s) poroelastic
integer:: Usolidnorm_index(1)
! maximum speed in velocity model
@@ -367,6 +378,13 @@
integer, dimension(:), allocatable :: coupling_ac_el_ispec
integer :: num_coupling_ac_el_faces
+! acoustic-poroelastic coupling surface
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: coupling_ac_po_normal
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: coupling_ac_po_jacobian2Dw
+ integer, dimension(:,:,:), allocatable :: coupling_ac_po_ijk
+ integer, dimension(:), allocatable :: coupling_ac_po_ispec
+ integer :: num_coupling_ac_po_faces
+
! material flag
logical, dimension(:), allocatable :: ispec_is_acoustic
integer, dimension(:,:), allocatable :: phase_ispec_inner_acoustic
@@ -411,8 +429,32 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_solid_poroelastic,&
rmass_fluid_poroelastic
+! displacement, velocity, acceleration
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accels_poroelastic,velocs_poroelastic,displs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accelw_poroelastic,velocw_poroelastic,displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accels_poroelastic,b_velocs_poroelastic,b_displs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accelw_poroelastic,b_velocw_poroelastic,b_displw_poroelastic
+
+! material properties
+! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: mustore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: etastore,tortstore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: phistore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: kappaarraystore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: permstore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vpI,rho_vpII,rho_vsI
+
+! 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 :: num_coupling_el_po_faces
+
! material flag
logical, dimension(:), allocatable :: ispec_is_poroelastic
+ integer, dimension(:,:), allocatable :: phase_ispec_inner_poroelastic
+ integer :: num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic
logical :: POROELASTIC_SIMULATION
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2011-10-06 02:26:04 UTC (rev 19024)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2011-10-06 02:32:28 UTC (rev 19025)
@@ -97,6 +97,17 @@
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
endif ! acoustic
+ ! poroelastic wave field
+ if( ispec_is_poroelastic(ispec) ) then
+ ! interpolates displ/veloc/accel at receiver locations
+ ! call compute_interpolated_dva(displw_poroelastic,velocw_poroelastic,accelw_poroelastic,NGLOB_AB, &
+ call compute_interpolated_dva(displs_poroelastic,velocs_poroelastic,accels_poroelastic,NGLOB_AB, &
+ ispec,NSPEC_AB,ibool, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ hxir,hetar,hgammar, &
+ dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
+ endif !poroelastic
+
!adjoint simulations
else if (SIMULATION_TYPE == 2) then
More information about the CIG-COMMITS
mailing list