[cig-commits] r20376 - in seismo/3D/SPECFEM3D_GLOBE/trunk/src: meshfem3D specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Fri Jun 15 05:56:43 PDT 2012
Author: joseph.charles
Date: 2012-06-15 05:56:42 -0700 (Fri, 15 Jun 2012)
New Revision: 20376
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
Log:
fixes Stacey boundary conditions bug by adding a C*deltat/2 contribution to the mass matrices on Stacey edges
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -29,65 +29,167 @@
nspec_actually,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore, &
gammaxstore,gammaystore,gammazstore, &
- iregion_code,nglob,rmass,rhostore,kappavstore, &
- nglob_oceans,rmass_ocean_load,NSPEC2D_TOP,ibelm_top,jacobian2D_top, &
- xstore,ystore,zstore,RHO_OCEANS)
+ iregion_code,rhostore,kappavstore, &
+ nglob_xy,nglob,prname, &
+ rmassx,rmassy,rmassz,DT,ichunk,NCHUNKS,ABSORBING_CONDITIONS, &
+ nglob_oceans,rmass_ocean_load, &
+ xstore,ystore,zstore,RHO_OCEANS, &
+ NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax, &
+ rho_vp,rho_vs,nspec_stacey, &
+ jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
+ jacobian2D_bottom,jacobian2D_top)
-! creates rmass and rmass_ocean_load
+ ! creates rmassx, rmassy, rmassz and rmass_ocean_load
use meshfem3D_models_par
implicit none
- integer myrank,nspec
+ integer :: myrank,nspec
+ integer :: idoubling(nspec)
+ integer :: ibool(NGLLX,NGLLY,NGLLZ,nspec)
+ integer :: nspec_actually
+ integer :: ichunk,NCHUNKS
- integer idoubling(nspec)
+ logical :: ABSORBING_CONDITIONS
- double precision wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ)
-
- integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
- integer nspec_actually
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_actually) :: &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore
- integer iregion_code
+ integer :: iregion_code,nglob_xy,nglob
- ! mass matrix
- integer nglob
- real(kind=CUSTOM_REAL), dimension(nglob) :: rmass
+ ! mass matrices
+ ! add C delta/2 contribution to the mass matrices on the Stacey edges
+ real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rhostore,kappavstore
! ocean mass matrix
- integer nglob_oceans
+ integer :: nglob_oceans
real(kind=CUSTOM_REAL), dimension(nglob_oceans) :: rmass_ocean_load
- integer NSPEC2D_TOP
- integer, dimension(NSPEC2D_TOP) :: ibelm_top
+ ! arrays with the mesh in double precision
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ystore
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: zstore
+
+ double precision :: RHO_OCEANS
+
+ ! processor identification
+ character(len=150) prname
+
+ ! time scheme
+ double precision :: DT,deltat,deltatover2
+
+ ! Stacey conditions put back
+ integer :: NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
+ integer :: nspec_stacey
+
+ double precision, dimension(NGLLX) :: wxgll
+ double precision, dimension(NGLLY) :: wygll
+ double precision, dimension(NGLLZ) :: wzgll
+
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX) :: jacobian2D_xmin,jacobian2D_xmax
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX) :: jacobian2D_ymin,jacobian2D_ymax
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_BOTTOM) :: jacobian2D_bottom
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_TOP) :: jacobian2D_top
- ! arrays with the mesh in double precision
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX) :: normal_xmin,normal_xmax
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX) :: normal_ymin,normal_ymax
- double precision RHO_OCEANS
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_stacey) :: rho_vp,rho_vs
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+ real(kind=CUSTOM_REAL) :: tx,ty,tz,sn
+ real(kind=CUSTOM_REAL) :: nx,ny,nz,vn
+
+ integer, dimension(NSPEC2D_TOP) :: ibelm_top
+ integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
+ integer, dimension(NSPEC2DMAX_XMIN_XMAX) :: ibelm_xmin,ibelm_xmax
+ integer, dimension(NSPEC2DMAX_YMIN_YMAX) :: ibelm_ymin,ibelm_ymax
+
+ integer, dimension(2,NSPEC2DMAX_YMIN_YMAX) :: nimin,nimax,nkmin_eta
+ integer, dimension(2,NSPEC2DMAX_XMIN_XMAX) :: njmin,njmax,nkmin_xi
+
+ integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,ispec2D
+
! local parameters
- double precision weight
- double precision xval,yval,zval,rval,thetaval,phival
- double precision lat,lon,colat
- double precision elevation,height_oceans
+ double precision :: xval,yval,zval,rval,thetaval,phival,weight
+ double precision :: lat,lon,colat
+ double precision :: elevation,height_oceans
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
- integer :: ispec,i,j,k,iglobnum
+ integer :: ispec,i,j,k,iglob
integer :: ix_oceans,iy_oceans,iz_oceans,ispec_oceans,ispec2D_top_crust
+ ! initializes matrices
+ !
+ ! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
- ! initializes
- rmass(:) = 0._CUSTOM_REAL
+ rmassx(:) = 0._CUSTOM_REAL
+ rmassy(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
+ ! use the non-dimensional time step to make the mass matrix correction
+ deltat = DT*dsqrt(PI*GRAV*RHOAV)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ do i=1,NGLLX
+ do j=1,NGLLY
+ wgllwgll_xy(i,j) = sngl(wxgll(i)*wygll(j))
+ enddo
+ enddo
+
+ do i=1,NGLLX
+ do k=1,NGLLZ
+ wgllwgll_xz(i,k) = sngl(wxgll(i)*wzgll(k))
+ enddo
+ enddo
+
+ do j=1,NGLLY
+ do k=1,NGLLZ
+ wgllwgll_yz(j,k) = sngl(wygll(j)*wzgll(k))
+ enddo
+ enddo
+
+ deltatover2 = sngl(0.5d0*deltat)
+
+ else ! double precision version
+ do i=1,NGLLX
+ do j=1,NGLLY
+ wgllwgll_xy(i,j) = wxgll(i)*wygll(j)
+ enddo
+ enddo
+
+ do i=1,NGLLX
+ do k=1,NGLLZ
+ wgllwgll_xz(i,k) = wxgll(i)*wzgll(k)
+ enddo
+ enddo
+
+ do j=1,NGLLY
+ do k=1,NGLLZ
+ wgllwgll_yz(j,k) = wygll(j)*wzgll(k)
+ enddo
+ enddo
+
+ deltatover2 = 0.5d0*deltat
+
+ endif
+
do ispec=1,nspec
! suppress fictitious elements in central cube
@@ -98,7 +200,7 @@
do i = 1,NGLLX
weight = wxgll(i)*wygll(j)*wzgll(k)
- iglobnum = ibool(i,j,k,ispec)
+ iglob = ibool(i,j,k,ispec)
! compute the jacobian
xixl = xixstore(i,j,k,ispec)
@@ -116,34 +218,36 @@
+ xizl*(etaxl*gammayl-etayl*gammaxl))
! definition depends if region is fluid or solid
- if(iregion_code == IREGION_CRUST_MANTLE .or. iregion_code == IREGION_INNER_CORE) then
+ select case( iregion_code)
+ case( IREGION_CRUST_MANTLE, IREGION_INNER_CORE )
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglobnum) = rmass(iglobnum) + &
+ rmassz(iglob) = rmassz(iglob) + &
sngl(dble(rhostore(i,j,k,ispec)) * dble(jacobianl) * weight)
else
- rmass(iglobnum) = rmass(iglobnum) + rhostore(i,j,k,ispec) * jacobianl * weight
+ rmassz(iglob) = rmassz(iglob) + rhostore(i,j,k,ispec) * jacobianl * weight
endif
! fluid in outer core
- else if(iregion_code == IREGION_OUTER_CORE) then
+ case( IREGION_OUTER_CORE )
! no anisotropy in the fluid, use kappav
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglobnum) = rmass(iglobnum) + &
+ rmassz(iglob) = rmassz(iglob) + &
sngl(dble(jacobianl) * weight * dble(rhostore(i,j,k,ispec)) / dble(kappavstore(i,j,k,ispec)))
else
- rmass(iglobnum) = rmass(iglobnum) + &
+ rmassz(iglob) = rmassz(iglob) + &
jacobianl * weight * rhostore(i,j,k,ispec) / kappavstore(i,j,k,ispec)
endif
- else
+ case default
call exit_MPI(myrank,'wrong region code')
- endif
+ end select
+
enddo
enddo
enddo
@@ -166,7 +270,7 @@
do ix_oceans = 1,NGLLX
do iy_oceans = 1,NGLLY
- iglobnum=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
+ iglob=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
! if 3D Earth, compute local height of oceans
if(CASE_3D) then
@@ -205,14 +309,15 @@
endif
! take into account inertia of water column
- weight = wxgll(ix_oceans)*wygll(iy_oceans)*dble(jacobian2D_top(ix_oceans,iy_oceans,ispec2D_top_crust)) &
- * dble(RHO_OCEANS) * height_oceans
+ weight = wxgll(ix_oceans) * wygll(iy_oceans) &
+ * dble(jacobian2D_top(ix_oceans,iy_oceans,ispec2D_top_crust)) &
+ * dble(RHO_OCEANS) * height_oceans
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + sngl(weight)
+ rmass_ocean_load(iglob) = rmass_ocean_load(iglob) + sngl(weight)
else
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + weight
+ rmass_ocean_load(iglob) = rmass_ocean_load(iglob) + weight
endif
enddo
@@ -221,8 +326,355 @@
enddo
! add regular mass matrix to ocean load contribution
- rmass_ocean_load(:) = rmass_ocean_load(:) + rmass(:)
+ rmass_ocean_load(:) = rmass_ocean_load(:) + rmassz(:)
endif
+ ! add C*delta/2 contribution to the mass matrices on the Stacey edges
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+
+ ! read arrays for Stacey conditions
+ open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
+ status='old',form='unformatted',action='read')
+ read(27) nimin
+ read(27) nimax
+ read(27) njmin
+ read(27) njmax
+ read(27) nkmin_xi
+ read(27) nkmin_eta
+ close(27)
+
+ select case(iregion_code)
+
+ case(IREGION_CRUST_MANTLE)
+
+ rmassx(:) = rmassz(:)
+ rmassy(:) = rmassz(:)
+
+ ! xmin
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AC) then
+
+ do ispec2D=1,nspec2D_xmin
+
+ ispec=ibelm_xmin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_xmin(1,j,k,ispec2D)
+ ny = normal_xmin(2,j,k,ispec2D)
+ nz = normal_xmin(3,j,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AC
+
+ ! xmax
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AB) then
+
+ do ispec2D=1,nspec2D_xmax
+
+ ispec=ibelm_xmax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_xmax(1,j,k,ispec2D)
+ ny = normal_xmax(2,j,k,ispec2D)
+ nz = normal_xmax(3,j,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AB
+
+ ! ymin
+ do ispec2D=1,nspec2D_ymin
+
+ ispec=ibelm_ymin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_ymin(1,i,k,ispec2D)
+ ny = normal_ymin(2,i,k,ispec2D)
+ nz = normal_ymin(3,i,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! ymax
+ do ispec2D=1,nspec2D_ymax
+
+ ispec=ibelm_ymax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ nx = normal_ymax(1,i,k,ispec2D)
+ ny = normal_ymax(2,i,k,ispec2D)
+ nz = normal_ymax(3,i,k,ispec2D)
+
+ vn = deltatover2*(nx+ny+nz)
+
+ tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
+ ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
+ tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
+
+ weight = jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
+ rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
+ rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ else
+ rmassx(iglob) = rmassx(iglob) + tx*weight
+ rmassy(iglob) = rmassy(iglob) + ty*weight
+ rmassz(iglob) = rmassz(iglob) + tz*weight
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! check that mass matrix is positive
+ if(minval(rmassx(:)) <= 0.) call exit_MPI(myrank,'negative rmassx matrix term')
+ if(minval(rmassy(:)) <= 0.) call exit_MPI(myrank,'negative rmassy matrix term')
+
+ case(IREGION_OUTER_CORE)
+
+ ! xmin
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AC) then
+
+ do ispec2D=1,nspec2D_xmin
+
+ ispec=ibelm_xmin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AC
+
+ ! xmax
+ ! if two chunks exclude this face for one of them
+ if(NCHUNKS == 1 .or. ichunk == CHUNK_AB) then
+
+ do ispec2D=1,nspec2D_xmax
+
+ ispec=ibelm_xmax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ endif ! NCHUNKS == 1 .or. ichunk == CHUNK_AB
+
+ ! ymin
+ do ispec2D=1,nspec2D_ymin
+
+ ispec=ibelm_ymin(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! ymax
+ do ispec2D=1,nspec2D_ymax
+
+ ispec=ibelm_ymax(ispec2D)
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! bottom (zmin)
+ do ispec2D=1,NSPEC2D_BOTTOM
+
+ ispec=ibelm_bottom(ispec2D)
+
+ k=1
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob=ibool(i,j,k,ispec)
+
+ sn = deltatover2/rho_vp(i,j,k,ispec)
+
+ weight = jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassz(iglob) = rmassz(iglob) + sngl(weight*sn)
+ else
+ rmassz(iglob) = rmassz(iglob) + weight*sn
+ endif
+ enddo
+ enddo
+ enddo
+
+ case( IREGION_INNER_CORE )
+
+ case default
+ call exit_MPI(myrank,'wrong region code')
+
+ end select
+
+ endif
+
+ ! check that mass matrix is positive
+ if(minval(rmassz(:)) <= 0.) call exit_MPI(myrank,'negative rmassz matrix term')
+
end subroutine create_mass_matrices
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -165,8 +165,9 @@
integer nglob,nglob_theor,ieoff,ilocnum,ier
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
+ ! mass matrices
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
+ integer :: nglob_xy
! mass matrix and bathymetry for ocean load
integer nglob_oceans
@@ -997,9 +998,33 @@
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,&
xigll,yigll,zigll)
- ! allocates mass matrix in this slice (will be fully assembled in the solver)
- allocate(rmass(nglob),stat=ier)
+ ! allocates mass matrices in this slice (will be fully assembled in the solver)
+ !
+ ! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+ select case(iregion_code)
+ case( IREGION_CRUST_MANTLE )
+ nglob_xy = nglob
+ case( IREGION_INNER_CORE, IREGION_OUTER_CORE )
+ nglob_xy = 1
+ endselect
+ else
+ nglob_xy = 1
+ endif
+
+ allocate(rmassx(nglob_xy),stat=ier)
if(ier /= 0) stop 'error in allocate 21'
+ allocate(rmassy(nglob_xy),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+ allocate(rmassz(nglob),stat=ier)
+ if(ier /= 0) stop 'error in allocate 21'
+
! allocates ocean load mass matrix as well if oceans
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) then
nglob_oceans = nglob
@@ -1010,14 +1035,23 @@
allocate(rmass_ocean_load(nglob_oceans),stat=ier)
if(ier /= 0) stop 'error in allocate 22'
- ! creating mass matrix in this slice (will be fully assembled in the solver)
+ ! creating mass matrices in this slice (will be fully assembled in the solver)
call create_mass_matrices(myrank,nspec,idoubling,wxgll,wygll,wzgll,ibool, &
nspec_actually,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore, &
gammaxstore,gammaystore,gammazstore, &
- iregion_code,nglob,rmass,rhostore,kappavstore, &
- nglob_oceans,rmass_ocean_load,NSPEC2D_TOP,ibelm_top,jacobian2D_top, &
- xstore,ystore,zstore,RHO_OCEANS)
+ iregion_code,rhostore,kappavstore, &
+ nglob_xy,nglob,prname, &
+ rmassx,rmassy,rmassz,DT,ichunk,NCHUNKS,ABSORBING_CONDITIONS, &
+ nglob_oceans,rmass_ocean_load, &
+ xstore,ystore,zstore,RHO_OCEANS, &
+ NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax, &
+ rho_vp,rho_vs,nspec_stacey, &
+ jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
+ jacobian2D_bottom,jacobian2D_top)
! save the binary files
#ifdef USE_SERIAL_CASCADE_FOR_IOs
@@ -1034,25 +1068,28 @@
nspec_ani,c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,is_on_a_slice_edge,rmass,rmass_ocean_load,nglob_oceans, &
+ ibool,idoubling,is_on_a_slice_edge,nglob_xy,nglob, &
+ rmassx,rmassy,rmassz,rmass_ocean_load,nglob_oceans, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top, &
jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top,nspec,nglob, &
+ jacobian2D_bottom,jacobian2D_top,nspec, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION, &
size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4),size(tau_e_store,5),&
- ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
+ NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
#ifdef USE_SERIAL_CASCADE_FOR_IOs
you_can_start_doing_IOs = .true.
if (myrank < NPROC_XI*NPROC_ETA-1) call MPI_SEND(you_can_start_doing_IOs, 1, MPI_LOGICAL, myrank+1, itag, MPI_COMM_WORLD, ier)
#endif
- deallocate(rmass,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(rmass_ocean_load,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(rmassx,stat=ier); if(ier /= 0) stop 'error in deallocate rmassx'
+ deallocate(rmassy,stat=ier); if(ier /= 0) stop 'error in deallocate rmassy'
+ deallocate(rmassz,stat=ier); if(ier /= 0) stop 'error in deallocate rmassz'
+ deallocate(rmass_ocean_load,stat=ier); if(ier /= 0) stop 'error in deallocate rmass_ocean_load'
! boundary mesh
if (SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -33,19 +33,19 @@
nspec_ani,c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,is_on_a_slice_edge,rmass_and_also_temporary_array,rmass_ocean_load,npointot_oceans, &
+ ibool,idoubling,is_on_a_slice_edge,nglob_xy,nglob, &
+ rmassx,rmassy,rmassz,rmass_ocean_load,npointot_oceans, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top, &
jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top,nspec,nglob, &
+ jacobian2D_bottom,jacobian2D_top,nspec, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION,vx,vy,vz,vnspec, &
- ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
+ NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
-
implicit none
include "constants.h"
@@ -53,9 +53,9 @@
logical ATTENUATION
character(len=150) prname
- integer iregion_code
+ integer iregion_code,NCHUNKS
- integer nspec,nglob,nspec_stacey
+ integer nspec,nglob_xy,nglob,nspec_stacey
integer npointot_oceans
! Stacey
@@ -91,11 +91,11 @@
! this for non blocking MPI
logical, dimension(nspec) :: is_on_a_slice_edge
-! mass matrix
-! use rmass_and_also_temporary_array as temporary storage for the rest of this routine
-! once it is written to disk and thus not needed any more;
-! thus use a variable name that underlines this
- real(kind=CUSTOM_REAL) rmass_and_also_temporary_array(nglob)
+! mass matrices
+! use rmassz as temporary storage for the rest of this routine
+! once it is written to disk and thus not needed any more
+ real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
! additional ocean load mass matrix
real(kind=CUSTOM_REAL) rmass_ocean_load(npointot_oceans)
@@ -247,8 +247,20 @@
endif
-! mass matrix
- write(27) rmass_and_also_temporary_array
+! mass matrices
+!
+! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+!
+! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) then
+ write(27) rmassx
+ write(27) rmassy
+ endif
+
+ write(27) rmassz
! additional ocean load mass matrix if oceans and if we are in the crust
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) write(27) rmass_ocean_load
@@ -256,12 +268,13 @@
close(27)
open(unit=27,file=prname(1:len_trim(prname))//'solver_data_2.bin',status='unknown',form='unformatted',action='write')
+
! mesh arrays used in the solver to locate source and receivers
! and for anisotropy and gravity, save in single precision
-! use rmass_and_also_temporary_array for temporary storage to perform conversion, since already saved
+! use rmassz for temporary storage to perform conversion, since already saved
!--- x coordinate
- rmass_and_also_temporary_array(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -269,18 +282,18 @@
iglob = ibool(i,j,k,ispec)
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass_and_also_temporary_array(iglob) = sngl(xstore(i,j,k,ispec))
+ rmassz(iglob) = sngl(xstore(i,j,k,ispec))
else
- rmass_and_also_temporary_array(iglob) = xstore(i,j,k,ispec)
+ rmassz(iglob) = xstore(i,j,k,ispec)
endif
enddo
enddo
enddo
enddo
- write(27) rmass_and_also_temporary_array
+ write(27) rmassz
!--- y coordinate
- rmass_and_also_temporary_array(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -288,18 +301,18 @@
iglob = ibool(i,j,k,ispec)
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass_and_also_temporary_array(iglob) = sngl(ystore(i,j,k,ispec))
+ rmassz(iglob) = sngl(ystore(i,j,k,ispec))
else
- rmass_and_also_temporary_array(iglob) = ystore(i,j,k,ispec)
+ rmassz(iglob) = ystore(i,j,k,ispec)
endif
enddo
enddo
enddo
enddo
- write(27) rmass_and_also_temporary_array
+ write(27) rmassz
!--- z coordinate
- rmass_and_also_temporary_array(:) = 0._CUSTOM_REAL
+ rmassz(:) = 0._CUSTOM_REAL
do ispec = 1,nspec
do k = 1,NGLLZ
do j = 1,NGLLY
@@ -307,15 +320,15 @@
iglob = ibool(i,j,k,ispec)
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- rmass_and_also_temporary_array(iglob) = sngl(zstore(i,j,k,ispec))
+ rmassz(iglob) = sngl(zstore(i,j,k,ispec))
else
- rmass_and_also_temporary_array(iglob) = zstore(i,j,k,ispec)
+ rmassz(iglob) = zstore(i,j,k,ispec)
endif
enddo
enddo
enddo
enddo
- write(27) rmass_and_also_temporary_array
+ write(27) rmassz
write(27) ibool
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -439,22 +439,35 @@
!
subroutine compute_coupling_ocean(accel_crust_mantle,b_accel_crust_mantle, &
- rmass_crust_mantle,rmass_ocean_load,normal_top_crust_mantle, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load, &
- SIMULATION_TYPE,nspec_top)
+ updated_dof_ocean_load,NGLOB_XY, &
+ SIMULATION_TYPE,nspec_top, &
+ ABSORBING_CONDITIONS)
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
- accel_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE_ADJOINT) :: &
- b_accel_crust_mantle
+ integer :: NGLOB_XY
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: accel_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE_ADJOINT) :: b_accel_crust_mantle
+
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
+
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_CM) :: normal_top_crust_mantle
@@ -462,12 +475,15 @@
integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
logical, dimension(NGLOB_CRUST_MANTLE_OCEANS) :: updated_dof_ocean_load
+ logical :: ABSORBING_CONDITIONS
integer SIMULATION_TYPE
integer nspec_top
! local parameters
real(kind=CUSTOM_REAL) :: force_normal_comp,b_force_normal_comp
+ real(kind=CUSTOM_REAL) :: additional_term_x,additional_term_y,additional_term_z
+ real(kind=CUSTOM_REAL) :: b_additional_term_x,b_additional_term_y,b_additional_term_z
real(kind=CUSTOM_REAL) :: additional_term,b_additional_term
real(kind=CUSTOM_REAL) :: nx,ny,nz
integer :: i,j,k,ispec,ispec2D,iglob
@@ -475,61 +491,128 @@
! initialize the updates
updated_dof_ocean_load(:) = .false.
- ! for surface elements exactly at the top of the crust (ocean bottom)
- do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
- ispec = ibelm_top_crust_mantle(ispec2D)
+ ! for surface elements exactly at the top of the crust (ocean bottom)
+ do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
- ! only for DOFs exactly at the top of the crust (ocean bottom)
- k = NGLLZ
+ ispec = ibelm_top_crust_mantle(ispec2D)
- do j = 1,NGLLY
- do i = 1,NGLLX
+ ! only for DOFs exactly at the top of the crust (ocean bottom)
+ k = NGLLZ
- ! get global point number
- iglob = ibool_crust_mantle(i,j,k,ispec)
+ do j = 1,NGLLY
+ do i = 1,NGLLX
- ! only update once
- if(.not. updated_dof_ocean_load(iglob)) then
+ ! get global point number
+ iglob = ibool_crust_mantle(i,j,k,ispec)
- ! get normal
- nx = normal_top_crust_mantle(1,i,j,ispec2D)
- ny = normal_top_crust_mantle(2,i,j,ispec2D)
- nz = normal_top_crust_mantle(3,i,j,ispec2D)
+ ! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
- ! make updated component of right-hand side
- ! we divide by rmass_crust_mantle() which is 1 / M
- ! we use the total force which includes the Coriolis term above
- force_normal_comp = (accel_crust_mantle(1,iglob)*nx + &
- accel_crust_mantle(2,iglob)*ny + &
- accel_crust_mantle(3,iglob)*nz) / rmass_crust_mantle(iglob)
+ ! get normal
+ nx = normal_top_crust_mantle(1,i,j,ispec2D)
+ ny = normal_top_crust_mantle(2,i,j,ispec2D)
+ nz = normal_top_crust_mantle(3,i,j,ispec2D)
- additional_term = (rmass_ocean_load(iglob) - rmass_crust_mantle(iglob)) * force_normal_comp
+ ! make updated component of right-hand side
+ ! we divide by rmass_crust_mantle() which is 1 / M
+ ! we use the total force which includes the Coriolis term above
- accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term * nx
- accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term * ny
- accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term * nz
+ force_normal_comp = accel_crust_mantle(1,iglob)*nx / rmassx_crust_mantle(iglob) + &
+ accel_crust_mantle(2,iglob)*ny / rmassy_crust_mantle(iglob) + &
+ accel_crust_mantle(3,iglob)*nz / rmassz_crust_mantle(iglob)
- if (SIMULATION_TYPE == 3) then
- b_force_normal_comp = (b_accel_crust_mantle(1,iglob)*nx + &
- b_accel_crust_mantle(2,iglob)*ny + &
- b_accel_crust_mantle(3,iglob)*nz) / rmass_crust_mantle(iglob)
+ additional_term_x = (rmass_ocean_load(iglob) - rmassx_crust_mantle(iglob)) * force_normal_comp
+ additional_term_y = (rmass_ocean_load(iglob) - rmassy_crust_mantle(iglob)) * force_normal_comp
+ additional_term_z = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
- b_additional_term = (rmass_ocean_load(iglob) - rmass_crust_mantle(iglob)) * b_force_normal_comp
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term_x * nx
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term_y * ny
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term_z * nz
- b_accel_crust_mantle(1,iglob) = b_accel_crust_mantle(1,iglob) + b_additional_term * nx
- b_accel_crust_mantle(2,iglob) = b_accel_crust_mantle(2,iglob) + b_additional_term * ny
- b_accel_crust_mantle(3,iglob) = b_accel_crust_mantle(3,iglob) + b_additional_term * nz
- endif
+ if (SIMULATION_TYPE == 3) then
+ b_force_normal_comp = b_accel_crust_mantle(1,iglob)*nx / rmassx_crust_mantle(iglob) + &
+ b_accel_crust_mantle(2,iglob)*ny / rmassy_crust_mantle(iglob) + &
+ b_accel_crust_mantle(3,iglob)*nz / rmassz_crust_mantle(iglob)
- ! done with this point
- updated_dof_ocean_load(iglob) = .true.
+ b_additional_term_x = (rmass_ocean_load(iglob) - rmassx_crust_mantle(iglob)) * b_force_normal_comp
+ b_additional_term_y = (rmass_ocean_load(iglob) - rmassy_crust_mantle(iglob)) * b_force_normal_comp
+ b_additional_term_z = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * b_force_normal_comp
- endif
+ b_accel_crust_mantle(1,iglob) = b_accel_crust_mantle(1,iglob) + b_additional_term_x * nx
+ b_accel_crust_mantle(2,iglob) = b_accel_crust_mantle(2,iglob) + b_additional_term_y * ny
+ b_accel_crust_mantle(3,iglob) = b_accel_crust_mantle(3,iglob) + b_additional_term_z * nz
+ endif
- enddo
- enddo
- enddo
+ ! done with this point
+ updated_dof_ocean_load(iglob) = .true.
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ else
+
+ ! for surface elements exactly at the top of the crust (ocean bottom)
+ do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+
+ ispec = ibelm_top_crust_mantle(ispec2D)
+
+ ! only for DOFs exactly at the top of the crust (ocean bottom)
+ k = NGLLZ
+
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ ! get global point number
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+
+ ! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
+
+ ! get normal
+ nx = normal_top_crust_mantle(1,i,j,ispec2D)
+ ny = normal_top_crust_mantle(2,i,j,ispec2D)
+ nz = normal_top_crust_mantle(3,i,j,ispec2D)
+
+ ! make updated component of right-hand side
+ ! we divide by rmass_crust_mantle() which is 1 / M
+ ! we use the total force which includes the Coriolis term above
+ force_normal_comp = (accel_crust_mantle(1,iglob)*nx + &
+ accel_crust_mantle(2,iglob)*ny + &
+ accel_crust_mantle(3,iglob)*nz) / rmassz_crust_mantle(iglob)
+
+ additional_term = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
+
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term * nx
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term * ny
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term * nz
+
+ if (SIMULATION_TYPE == 3) then
+ b_force_normal_comp = (b_accel_crust_mantle(1,iglob)*nx + &
+ b_accel_crust_mantle(2,iglob)*ny + &
+ b_accel_crust_mantle(3,iglob)*nz) / rmassz_crust_mantle(iglob)
+
+ b_additional_term = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * b_force_normal_comp
+
+ b_accel_crust_mantle(1,iglob) = b_accel_crust_mantle(1,iglob) + b_additional_term * nx
+ b_accel_crust_mantle(2,iglob) = b_accel_crust_mantle(2,iglob) + b_additional_term * ny
+ b_accel_crust_mantle(3,iglob) = b_accel_crust_mantle(3,iglob) + b_additional_term * nz
+ endif
+
+ ! done with this point
+ updated_dof_ocean_load(iglob) = .true.
+
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ endif
+
end subroutine compute_coupling_ocean
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -25,8 +25,8 @@
!
!=====================================================================
- subroutine prepare_timerun_rmass(myrank,rmass_ocean_load,rmass_crust_mantle, &
- rmass_outer_core,rmass_inner_core, &
+ subroutine prepare_timerun_rmass(myrank,rmass_ocean_load,rmassx_crust_mantle,rmassy_crust_mantle, &
+ rmassz_crust_mantle,rmass_outer_core,rmass_inner_core, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle, &
iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
@@ -44,7 +44,7 @@
iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
buffer_send_faces_scalar,buffer_received_faces_scalar, &
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
- NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS,NGLOB_XY,ABSORBING_CONDITIONS, &
NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB2DMAX_XY,npoin2D_max_all_CM_IC)
implicit none
@@ -54,9 +54,14 @@
include "OUTPUT_FILES/values_from_mesher.h"
integer myrank,npoin2D_max_all_CM_IC
+ integer NGLOB_XY
+ logical ABSORBING_CONDITIONS
+
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: rmass_outer_core
real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: rmass_inner_core
@@ -134,7 +139,8 @@
endif
! crust and mantle
- call assemble_MPI_scalar_block(myrank,rmass_crust_mantle,NGLOB_CRUST_MANTLE, &
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ call assemble_MPI_scalar_block(myrank,rmassx_crust_mantle,NGLOB_CRUST_MANTLE, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
@@ -147,6 +153,33 @@
NPROC_XI_VAL,NPROC_ETA_VAL,NGLOB1D_RADIAL(IREGION_CRUST_MANTLE), &
NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_XY,NCHUNKS_VAL)
+ call assemble_MPI_scalar_block(myrank,rmassy_crust_mantle,NGLOB_CRUST_MANTLE, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
+ npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
+ iboolfaces_crust_mantle,iboolcorner_crust_mantle, &
+ iprocfrom_faces,iprocto_faces,imsg_type, &
+ iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+ buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_max_all_CM_IC, &
+ buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+ NPROC_XI_VAL,NPROC_ETA_VAL,NGLOB1D_RADIAL(IREGION_CRUST_MANTLE), &
+ NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_XY,NCHUNKS_VAL)
+ endif
+
+ call assemble_MPI_scalar_block(myrank,rmassz_crust_mantle,NGLOB_CRUST_MANTLE, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
+ npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
+ iboolfaces_crust_mantle,iboolcorner_crust_mantle, &
+ iprocfrom_faces,iprocto_faces,imsg_type, &
+ iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+ buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_max_all_CM_IC, &
+ buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+ NPROC_XI_VAL,NPROC_ETA_VAL,NGLOB1D_RADIAL(IREGION_CRUST_MANTLE), &
+ NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_XY,NCHUNKS_VAL)
+
! outer core
call assemble_MPI_scalar_block(myrank,rmass_outer_core,NGLOB_OUTER_CORE, &
iproc_xi,iproc_eta,ichunk,addressing, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -35,9 +35,9 @@
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,ispec_is_tiso, &
- is_on_a_slice_edge,rmass,rmass_ocean_load,nspec,nglob, &
- READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY, &
+ ibool,idoubling,ispec_is_tiso,nglob_xy,nglob, &
+ rmassx,rmassy,rmassz,rmass_ocean_load,nspec, &
+ is_on_a_slice_edge,READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY, &
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS,LOCAL_PATH,ABSORBING_CONDITIONS)
implicit none
@@ -54,7 +54,7 @@
character(len=150) LOCAL_PATH
- integer nspec,nglob
+ integer nspec,nglob,nglob_xy
integer nspec_iso,nspec_tiso,nspec_ani
@@ -83,9 +83,12 @@
real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,nspec)
real(kind=CUSTOM_REAL) rho_vs(NGLLX,NGLLY,NGLLZ,nspec)
-! mass matrix and additional ocean load mass matrix
- real(kind=CUSTOM_REAL), dimension(nglob) :: rmass,rmass_ocean_load
+! mass matrices and additional ocean load mass matrix
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
+
! global addressing
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
@@ -173,9 +176,21 @@
endif
-! mass matrix
- read(IIN) rmass
+! mass matrices
+!
+! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+!
+! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) then
+ read(IIN) rmassx
+ read(IIN) rmassy
+ endif
+ read(IIN) rmassz
+
! read additional ocean load mass matrix
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) read(IIN) rmass_ocean_load
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -42,7 +42,8 @@
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
! -- idoubling_crust_mantle
- is_on_a_slice_edge_crust_mantle,rmass_crust_mantle,rmass_ocean_load, &
+ is_on_a_slice_edge_crust_mantle,rmass_ocean_load, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
vp_outer_core,xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -59,7 +60,7 @@
c33store_inner_core,c44store_inner_core, &
ibool_inner_core,idoubling_inner_core,ispec_is_tiso_inner_core, &
is_on_a_slice_edge_inner_core,rmass_inner_core, &
- ABSORBING_CONDITIONS,LOCAL_PATH)
+ ABSORBING_CONDITIONS,LOCAL_PATH,NGLOB_XY)
implicit none
@@ -101,8 +102,20 @@
! integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling_crust_mantle
logical, dimension(NSPEC_CRUST_MANTLE) :: ispec_is_tiso_crust_mantle
- ! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ integer :: NGLOB_XY, NGLOB_DUMMY
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
+
! additional mass matrix for ocean load
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
@@ -145,6 +158,7 @@
!local parameters
logical READ_KAPPA_MU,READ_TISO
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,1) :: dummy_array
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: dummy_rmass
integer, dimension(NSPEC_CRUST_MANTLE) :: dummy_i
! this for non blocking MPI
@@ -179,6 +193,7 @@
READ_KAPPA_MU = .true.
READ_TISO = .true.
endif
+
call read_arrays_solver(IREGION_CRUST_MANTLE,myrank, &
rho_vp_crust_mantle,rho_vs_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
@@ -195,12 +210,9 @@
c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
- ibool_crust_mantle,dummy_i, &
- ! --idoubling_crust_mantle,
- ispec_is_tiso_crust_mantle, &
- is_on_a_slice_edge_crust_mantle,rmass_crust_mantle,rmass_ocean_load, &
- NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE, &
- READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
+ ibool_crust_mantle,dummy_i,ispec_is_tiso_crust_mantle,NGLOB_XY,NGLOB_CRUST_MANTLE, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle,rmass_ocean_load,NSPEC_CRUST_MANTLE, &
+ is_on_a_slice_edge_crust_mantle,READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS)
! synchronizes processes
@@ -214,6 +226,9 @@
nspec_tiso = 1
nspec_ani = 1
+ NGLOB_DUMMY = 1
+ allocate(dummy_rmass(NGLOB_DUMMY))
+
call read_arrays_solver(IREGION_OUTER_CORE,myrank, &
vp_outer_core,dummy_array, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
@@ -230,10 +245,9 @@
dummy_array,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
- ibool_outer_core,idoubling_outer_core,ispec_is_tiso_outer_core, &
- is_on_a_slice_edge_outer_core,rmass_outer_core,rmass_ocean_load, &
- NSPEC_OUTER_CORE,NGLOB_OUTER_CORE, &
- READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
+ ibool_outer_core,idoubling_outer_core,ispec_is_tiso_outer_core,NGLOB_XY,NGLOB_OUTER_CORE, &
+ dummy_rmass,dummy_rmass,rmass_outer_core,rmass_ocean_load,NSPEC_OUTER_CORE, &
+ is_on_a_slice_edge_outer_core,READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS)
! synchronizes processes
@@ -267,10 +281,9 @@
dummy_array,dummy_array,dummy_array, &
c44store_inner_core,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
- ibool_inner_core,idoubling_inner_core,ispec_is_tiso_inner_core, &
- is_on_a_slice_edge_inner_core,rmass_inner_core,rmass_ocean_load, &
- NSPEC_INNER_CORE,NGLOB_INNER_CORE, &
- READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
+ ibool_inner_core,idoubling_inner_core,ispec_is_tiso_inner_core,NGLOB_XY,NGLOB_INNER_CORE, &
+ dummy_rmass,dummy_rmass,rmass_inner_core,rmass_ocean_load,NSPEC_INNER_CORE, &
+ is_on_a_slice_edge_inner_core,READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS)
! check that the number of points in this slice is correct
@@ -285,6 +298,8 @@
if(minval(ibool_inner_core(:,:,:,:)) /= 1 .or. maxval(ibool_inner_core(:,:,:,:)) /= NGLOB_INNER_CORE) &
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in inner core')
+ deallocate(dummy_rmass)
+
! synchronizes processes
call sync_all()
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2012-06-15 11:07:03 UTC (rev 20375)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2012-06-15 12:56:42 UTC (rev 20376)
@@ -537,9 +537,20 @@
! integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling_crust_mantle
logical, dimension(NSPEC_CRUST_MANTLE) :: ispec_is_tiso_crust_mantle
-! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+! mass matrices
+!
+! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+!
+! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
+ integer :: NGLOB_XY
+
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle
@@ -624,10 +635,10 @@
alpha_kl_outer_core
! approximate hessian
- real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: hess_kl_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: hess_kl_crust_mantle
! check for deviatoric kernel for outer core region
- real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: beta_kl_outer_core
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: beta_kl_outer_core
integer :: nspec_beta_kl_outer_core
logical,parameter:: deviatoric_outercore = .false.
@@ -1034,6 +1045,25 @@
!!!!!!! print *,'starting doing serialized I/Os on rank ',myrank
#endif
+ ! allocates mass matrices in this slice (will be fully assembled in the solver)
+ !
+ ! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ NGLOB_XY = NGLOB_CRUST_MANTLE
+ else
+ NGLOB_XY = 1
+ endif
+
+ allocate(rmassx_crust_mantle(NGLOB_XY),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating rmassx_crust_mantle')
+ allocate(rmassy_crust_mantle(NGLOB_XY),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating rmassy_crust_mantle')
+
call read_mesh_databases(myrank,rho_vp_crust_mantle,rho_vs_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
@@ -1051,7 +1081,8 @@
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
! -- idoubling_crust_mantle,
- is_on_a_slice_edge_crust_mantle,rmass_crust_mantle,rmass_ocean_load, &
+ is_on_a_slice_edge_crust_mantle,rmass_ocean_load, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
vp_outer_core,xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
@@ -1068,7 +1099,7 @@
c33store_inner_core,c44store_inner_core, &
ibool_inner_core,idoubling_inner_core,ispec_is_tiso_inner_core, &
is_on_a_slice_edge_inner_core,rmass_inner_core, &
- ABSORBING_CONDITIONS,LOCAL_PATH)
+ ABSORBING_CONDITIONS,LOCAL_PATH,NGLOB_XY)
! read 2-D addressing for summation between slices with MPI
call read_mesh_databases_addressing(myrank, &
@@ -1513,8 +1544,8 @@
endif
! the mass matrix needs to be assembled with MPI here once and for all
- call prepare_timerun_rmass(myrank,rmass_ocean_load,rmass_crust_mantle, &
- rmass_outer_core,rmass_inner_core, &
+ call prepare_timerun_rmass(myrank,rmass_ocean_load,rmassx_crust_mantle,rmassy_crust_mantle, &
+ rmassz_crust_mantle,rmass_outer_core,rmass_inner_core, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle, &
iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
@@ -1532,7 +1563,7 @@
iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
buffer_send_faces,buffer_received_faces, &
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
- NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS,NGLOB_XY,ABSORBING_CONDITIONS, &
NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB2DMAX_XY,npoin2D_max_all_CM_IC)
! mass matrix including central cube
@@ -1610,17 +1641,35 @@
if(OCEANS_VAL) then
if(minval(rmass_ocean_load) <= 0.) call exit_MPI(myrank,'negative mass matrix term for the oceans')
endif
- if(minval(rmass_crust_mantle) <= 0.) call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
- if(minval(rmass_inner_core) <= 0.) call exit_MPI(myrank,'negative mass matrix term for the inner core')
- if(minval(rmass_outer_core) <= 0.) call exit_MPI(myrank,'negative mass matrix term for the outer core')
+ ! add C*delta/2 contribution to the mass matrices on the Stacey edges
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if(minval(rmassx_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
+ if(minval(rmassy_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
+ endif
+
+ if(minval(rmassz_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the crust_mantle')
+ if(minval(rmass_inner_core) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the inner core')
+ if(minval(rmass_outer_core) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the outer core')
+
! for efficiency, invert final mass matrix once and for all on each slice
if(OCEANS_VAL) rmass_ocean_load = 1._CUSTOM_REAL / rmass_ocean_load
- rmass_crust_mantle = 1._CUSTOM_REAL / rmass_crust_mantle
+
+ ! add C*delta/2 contribution to the mass matrices on the Stacey edges
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ rmassx_crust_mantle = 1._CUSTOM_REAL / rmassx_crust_mantle
+ rmassy_crust_mantle = 1._CUSTOM_REAL / rmassy_crust_mantle
+ endif
+
+ rmassz_crust_mantle = 1._CUSTOM_REAL / rmassz_crust_mantle
rmass_outer_core = 1._CUSTOM_REAL / rmass_outer_core
rmass_inner_core = 1._CUSTOM_REAL / rmass_inner_core
-
! change x, y, z to r, theta and phi once and for all
! IMPROVE dangerous: old name kept (xstore ystore zstore) for new values
@@ -3400,53 +3449,106 @@
endif
endif ! end of assembling forces with the central cube
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+
#ifdef _HANDOPT
-! way 2:
- if(imodulo_NGLOB_CRUST_MANTLE4 >= 1) then
- do i=1,imodulo_NGLOB_CRUST_MANTLE4
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ ! way 2:
+ if(imodulo_NGLOB_CRUST_MANTLE4 >= 1) then
+ do i=1,imodulo_NGLOB_CRUST_MANTLE4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
+ do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB_CRUST_MANTLE,4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
- enddo
- endif
- do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB_CRUST_MANTLE,4
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+
+ accel_crust_mantle(1,i+1) = accel_crust_mantle(1,i+1)*rmassx_crust_mantle(i+1) &
+ + two_omega_earth*veloc_crust_mantle(2,i+1)
+ accel_crust_mantle(2,i+1) = accel_crust_mantle(2,i+1)*rmassy_crust_mantle(i+1) &
+ - two_omega_earth*veloc_crust_mantle(1,i+1)
+ accel_crust_mantle(3,i+1) = accel_crust_mantle(3,i+1)*rmassz_crust_mantle(i+1)
+
+ accel_crust_mantle(1,i+2) = accel_crust_mantle(1,i+2)*rmassx_crust_mantle(i+2) &
+ + two_omega_earth*veloc_crust_mantle(2,i+2)
+ accel_crust_mantle(2,i+2) = accel_crust_mantle(2,i+2)*rmassy_crust_mantle(i+2) &
+ - two_omega_earth*veloc_crust_mantle(1,i+2)
+ accel_crust_mantle(3,i+2) = accel_crust_mantle(3,i+2)*rmassz_crust_mantle(i+2)
+
+ accel_crust_mantle(1,i+3) = accel_crust_mantle(1,i+3)*rmassx_crust_mantle(i+3) &
+ + two_omega_earth*veloc_crust_mantle(2,i+3)
+ accel_crust_mantle(2,i+3) = accel_crust_mantle(2,i+3)*rmassy_crust_mantle(i+3) &
+ - two_omega_earth*veloc_crust_mantle(1,i+3)
+ accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmassz_crust_mantle(i+3)
+ enddo
+#else
+ ! way 1:
+ do i=1,NGLOB_CRUST_MANTLE
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+#endif
- accel_crust_mantle(1,i+1) = accel_crust_mantle(1,i+1)*rmass_crust_mantle(i+1) &
+ else
+
+#ifdef _HANDOPT
+ ! way 2:
+ if(imodulo_NGLOB_CRUST_MANTLE4 >= 1) then
+ do i=1,imodulo_NGLOB_CRUST_MANTLE4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
+ do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB,4
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + two_omega_earth*veloc_crust_mantle(2,i)
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - two_omega_earth*veloc_crust_mantle(1,i)
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+
+ accel_crust_mantle(1,i+1) = accel_crust_mantle(1,i+1)*rmassz_crust_mantle(i+1) &
+ two_omega_earth*veloc_crust_mantle(2,i+1)
- accel_crust_mantle(2,i+1) = accel_crust_mantle(2,i+1)*rmass_crust_mantle(i+1) &
+ accel_crust_mantle(2,i+1) = accel_crust_mantle(2,i+1)*rmassz_crust_mantle(i+1) &
- two_omega_earth*veloc_crust_mantle(1,i+1)
- accel_crust_mantle(3,i+1) = accel_crust_mantle(3,i+1)*rmass_crust_mantle(i+1)
+ accel_crust_mantle(3,i+1) = accel_crust_mantle(3,i+1)*rmassz_crust_mantle(i+1)
- accel_crust_mantle(1,i+2) = accel_crust_mantle(1,i+2)*rmass_crust_mantle(i+2) &
+ accel_crust_mantle(1,i+2) = accel_crust_mantle(1,i+2)*rmassz_crust_mantle(i+2) &
+ two_omega_earth*veloc_crust_mantle(2,i+2)
- accel_crust_mantle(2,i+2) = accel_crust_mantle(2,i+2)*rmass_crust_mantle(i+2) &
+ accel_crust_mantle(2,i+2) = accel_crust_mantle(2,i+2)*rmassz_crust_mantle(i+2) &
- two_omega_earth*veloc_crust_mantle(1,i+2)
- accel_crust_mantle(3,i+2) = accel_crust_mantle(3,i+2)*rmass_crust_mantle(i+2)
+ accel_crust_mantle(3,i+2) = accel_crust_mantle(3,i+2)*rmassz_crust_mantle(i+2)
- accel_crust_mantle(1,i+3) = accel_crust_mantle(1,i+3)*rmass_crust_mantle(i+3) &
+ accel_crust_mantle(1,i+3) = accel_crust_mantle(1,i+3)*rmassz_crust_mantle(i+3) &
+ two_omega_earth*veloc_crust_mantle(2,i+3)
- accel_crust_mantle(2,i+3) = accel_crust_mantle(2,i+3)*rmass_crust_mantle(i+3) &
+ accel_crust_mantle(2,i+3) = accel_crust_mantle(2,i+3)*rmassz_crust_mantle(i+3) &
- two_omega_earth*veloc_crust_mantle(1,i+3)
- accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmass_crust_mantle(i+3)
- enddo
+ accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmassz_crust_mantle(i+3)
+ enddo
#else
-! way 1:
- do i=1,NGLOB_CRUST_MANTLE
- accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ ! way 1:
+ do i=1,NGLOB_CRUST_MANTLE
+ accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ two_omega_earth*veloc_crust_mantle(2,i)
- accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ accel_crust_mantle(2,i) = accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
- two_omega_earth*veloc_crust_mantle(1,i)
- accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
- enddo
+ accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
#endif
+ endif
+
if (SIMULATION_TYPE == 3) then
! ------------------- new non blocking implementation -------------------
@@ -3715,64 +3817,118 @@
! ------------------- new non blocking implementation -------------------
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+
#ifdef _HANDOPT
-! way 2:
- if( imodulo_NGLOB_CRUST_MANTLE4 >=1 ) then
- do i=1,imodulo_NGLOB_CRUST_MANTLE4
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ ! way 2:
+ if( imodulo_NGLOB_CRUST_MANTLE4 >=1 ) then
+ do i=1,imodulo_NGLOB_CRUST_MANTLE4
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i)
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i)
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
+ do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB_CRUST_MANTLE,4
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
- b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
- b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmass_crust_mantle(i)
- enddo
- endif
- do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB_CRUST_MANTLE,4
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+
+ b_accel_crust_mantle(1,i+1) = b_accel_crust_mantle(1,i+1)*rmassx_crust_mantle(i+1) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i+1)
+ b_accel_crust_mantle(2,i+1) = b_accel_crust_mantle(2,i+1)*rmassy_crust_mantle(i+1) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i+1)
+ b_accel_crust_mantle(3,i+1) = b_accel_crust_mantle(3,i+1)*rmassz_crust_mantle(i+1)
+
+ b_accel_crust_mantle(1,i+2) = b_accel_crust_mantle(1,i+2)*rmassx_crust_mantle(i+2) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i+2)
+ b_accel_crust_mantle(2,i+2) = b_accel_crust_mantle(2,i+2)*rmassy_crust_mantle(i+2) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i+2)
+ b_accel_crust_mantle(3,i+2) = b_accel_crust_mantle(3,i+2)*rmassz_crust_mantle(i+2)
+
+ b_accel_crust_mantle(1,i+3) = b_accel_crust_mantle(1,i+3)*rmassx_crust_mantle(i+3) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i+3)
+ b_accel_crust_mantle(2,i+3) = b_accel_crust_mantle(2,i+3)*rmassy_crust_mantle(i+3) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i+3)
+ b_accel_crust_mantle(3,i+3) = b_accel_crust_mantle(3,i+3)*rmassz_crust_mantle(i+3)
+ enddo
+#else
+ ! way 1:
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
- b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassy_crust_mantle(i) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
- b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmass_crust_mantle(i)
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+#endif
- b_accel_crust_mantle(1,i+1) = b_accel_crust_mantle(1,i+1)*rmass_crust_mantle(i+1) &
+ else
+
+#ifdef _HANDOPT
+ ! way 2:
+ if( imodulo_NGLOB_CRUST_MANTLE4 >=1 ) then
+ do i=1,imodulo_NGLOB_CRUST_MANTLE4
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i)
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i)
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
+ endif
+ do i=imodulo_NGLOB_CRUST_MANTLE4+1,NGLOB_CRUST_MANTLE,4
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ + b_two_omega_earth*b_veloc_crust_mantle(2,i)
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
+ - b_two_omega_earth*b_veloc_crust_mantle(1,i)
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+
+ b_accel_crust_mantle(1,i+1) = b_accel_crust_mantle(1,i+1)*rmassz_crust_mantle(i+1) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i+1)
- b_accel_crust_mantle(2,i+1) = b_accel_crust_mantle(2,i+1)*rmass_crust_mantle(i+1) &
+ b_accel_crust_mantle(2,i+1) = b_accel_crust_mantle(2,i+1)*rmassz_crust_mantle(i+1) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i+1)
- b_accel_crust_mantle(3,i+1) = b_accel_crust_mantle(3,i+1)*rmass_crust_mantle(i+1)
+ b_accel_crust_mantle(3,i+1) = b_accel_crust_mantle(3,i+1)*rmassz_crust_mantle(i+1)
- b_accel_crust_mantle(1,i+2) = b_accel_crust_mantle(1,i+2)*rmass_crust_mantle(i+2) &
+ b_accel_crust_mantle(1,i+2) = b_accel_crust_mantle(1,i+2)*rmassz_crust_mantle(i+2) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i+2)
- b_accel_crust_mantle(2,i+2) = b_accel_crust_mantle(2,i+2)*rmass_crust_mantle(i+2) &
+ b_accel_crust_mantle(2,i+2) = b_accel_crust_mantle(2,i+2)*rmassz_crust_mantle(i+2) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i+2)
- b_accel_crust_mantle(3,i+2) = b_accel_crust_mantle(3,i+2)*rmass_crust_mantle(i+2)
+ b_accel_crust_mantle(3,i+2) = b_accel_crust_mantle(3,i+2)*rmassz_crust_mantle(i+2)
- b_accel_crust_mantle(1,i+3) = b_accel_crust_mantle(1,i+3)*rmass_crust_mantle(i+3) &
+ b_accel_crust_mantle(1,i+3) = b_accel_crust_mantle(1,i+3)*rmassz_crust_mantle(i+3) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i+3)
- b_accel_crust_mantle(2,i+3) = b_accel_crust_mantle(2,i+3)*rmass_crust_mantle(i+3) &
+ b_accel_crust_mantle(2,i+3) = b_accel_crust_mantle(2,i+3)*rmassz_crust_mantle(i+3) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i+3)
- b_accel_crust_mantle(3,i+3) = b_accel_crust_mantle(3,i+3)*rmass_crust_mantle(i+3)
- enddo
+ b_accel_crust_mantle(3,i+3) = b_accel_crust_mantle(3,i+3)*rmassz_crust_mantle(i+3)
+ enddo
#else
-! way 1:
- do i=1,NGLOB_CRUST_MANTLE
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmass_crust_mantle(i) &
+ ! way 1:
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ b_two_omega_earth*b_veloc_crust_mantle(2,i)
- b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmass_crust_mantle(i) &
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*rmassz_crust_mantle(i) &
- b_two_omega_earth*b_veloc_crust_mantle(1,i)
- b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmass_crust_mantle(i)
- enddo
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*rmassz_crust_mantle(i)
+ enddo
#endif
- endif ! SIMULATION_TYPE == 3
+ endif
+ endif ! SIMULATION_TYPE == 3
+
! couples ocean with crust mantle
- if(OCEANS_VAL) &
- call compute_coupling_ocean(accel_crust_mantle,b_accel_crust_mantle, &
- rmass_crust_mantle,rmass_ocean_load,normal_top_crust_mantle, &
- ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load, &
- SIMULATION_TYPE,NSPEC2D_TOP(IREGION_CRUST_MANTLE))
+ if(OCEANS_VAL) &
+ call compute_coupling_ocean(accel_crust_mantle,b_accel_crust_mantle, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ rmass_ocean_load,normal_top_crust_mantle, &
+ ibool_crust_mantle,ibelm_top_crust_mantle, &
+ updated_dof_ocean_load,NGLOB_XY, &
+ SIMULATION_TYPE,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+ ABSORBING_CONDITIONS)
-
!
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
@@ -4646,6 +4802,10 @@
noise_surface_movie)
endif
+ ! mass matrices
+ deallocate(rmassx_crust_mantle)
+ deallocate(rmassy_crust_mantle)
+
! close the main output file
if(myrank == 0) then
write(IMAIN,*)
More information about the CIG-COMMITS
mailing list