[cig-commits] r22574 - in seismo/3D/SPECFEM3D_GLOBE/trunk/src: meshfem3D specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Fri Jul 12 10:09:21 PDT 2013
Author: xie.zhinan
Date: 2013-07-12 10:09:21 -0700 (Fri, 12 Jul 2013)
New Revision: 22574
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/meshfem3D.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/initialize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.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:
finish changes concerning flag EXACT_MASS_MATRIX_FOR_ROTATION
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -40,7 +40,9 @@
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)
+ jacobian2D_bottom,jacobian2D_top,&
+ SIMULATION_TYPE,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK,&
+ nglob_xy_backward,b_rmassx,b_rmassy)
! creates rmassx, rmassy, rmassz and rmass_ocean_load
@@ -48,7 +50,7 @@
implicit none
- integer :: myrank,nspec
+ integer :: myrank,nspec,SIMULATION_TYPE
integer :: idoubling(nspec)
integer :: ibool(NGLLX,NGLLY,NGLLZ,nspec)
integer :: nspec_actually
@@ -65,7 +67,14 @@
! 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) :: two_omega_earth,scale_t_inv
+ ! mass matrices for backward simulation when SIMULATION_TYPE =3 and ROTATION is .true.
+ integer :: nglob_xy_backward
+ logical :: EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
+ real(kind=CUSTOM_REAL), dimension(nglob_xy_backward) :: b_rmassx,b_rmassy
+ real(kind=CUSTOM_REAL) :: b_two_omega_earth
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rhostore,kappavstore
! ocean mass matrix
@@ -140,9 +149,34 @@
rmassy(:) = 0._CUSTOM_REAL
rmassz(:) = 0._CUSTOM_REAL
+ b_rmassx(:) = 0._CUSTOM_REAL
+ b_rmassy(:) = 0._CUSTOM_REAL
+
! use the non-dimensional time step to make the mass matrix correction
deltat = DT*dsqrt(PI*GRAV*RHOAV)
+ if(ROTATION .and. (.not. USE_LDDRK) .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ ! distinguish between single and double precision for reals
+ scale_t_inv = dsqrt(PI*GRAV*RHOAV)
+ two_omega_earth = 0._CUSTOM_REAL
+ b_two_omega_earth = 0._CUSTOM_REAL
+ if (SIMULATION_TYPE /= 3) then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ two_omega_earth = sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat)
+ else
+ two_omega_earth = 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat
+ endif
+ endif
+
+ if (SIMULATION_TYPE == 3) then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ b_two_omega_earth = - sngl(2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat)
+ else
+ b_two_omega_earth = - 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv) * deltat
+ endif
+ endif
+ endif
+
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
do i=1,NGLLX
@@ -227,6 +261,33 @@
rmassz(iglob) = rmassz(iglob) + rhostore(i,j,k,ispec) * jacobianl * weight
endif
+ if(ROTATION .and. (.not. USE_LDDRK) .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ if(SIMULATION_TYPE == 3)then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ b_rmassx(iglob) = rmassz(iglob) - b_two_omega_earth * 0.5_CUSTOM_REAL*sngl(dble(jacobianl) * weight)
+ b_rmassy(iglob) = rmassz(iglob) + b_two_omega_earth * 0.5_CUSTOM_REAL*sngl(dble(jacobianl) * weight)
+ else
+ b_rmassx(iglob) = rmassz(iglob) - b_two_omega_earth * 0.5_CUSTOM_REAL * jacobianl * weight
+ b_rmassy(iglob) = rmassz(iglob) + b_two_omega_earth * 0.5_CUSTOM_REAL * jacobianl * weight
+ endif
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassz(iglob) - two_omega_earth * 0.5_CUSTOM_REAL*sngl(dble(jacobianl) * weight)
+ rmassy(iglob) = rmassz(iglob) + two_omega_earth * 0.5_CUSTOM_REAL*sngl(dble(jacobianl) * weight)
+ else
+ rmassx(iglob) = rmassz(iglob) - two_omega_earth * 0.5_CUSTOM_REAL * jacobianl * weight
+ rmassy(iglob) = rmassz(iglob) + two_omega_earth * 0.5_CUSTOM_REAL * jacobianl * weight
+ endif
+ else
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmassx(iglob) = rmassz(iglob) - two_omega_earth * 0.5_CUSTOM_REAL*sngl(dble(jacobianl) * weight)
+ rmassy(iglob) = rmassz(iglob) + two_omega_earth * 0.5_CUSTOM_REAL*sngl(dble(jacobianl) * weight)
+ else
+ rmassx(iglob) = rmassz(iglob) - two_omega_earth * 0.5_CUSTOM_REAL * jacobianl * weight
+ rmassy(iglob) = rmassz(iglob) + two_omega_earth * 0.5_CUSTOM_REAL * jacobianl * weight
+ endif
+ endif
+ endif
+
! fluid in outer core
case( IREGION_OUTER_CORE )
@@ -334,7 +395,7 @@
endif
! add C*deltat/2 contribution to the mass matrices on the Stacey edges
- if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS .and. (.not. USE_LDDRK)) then
! read arrays for Stacey conditions
open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
@@ -350,10 +411,12 @@
select case(iregion_code)
case(IREGION_CRUST_MANTLE)
+
+ if(.not. ROTATION)then
+ rmassx(:) = rmassz(:)
+ rmassy(:) = rmassz(:)
+ endif
- rmassx(:) = rmassz(:)
- rmassy(:) = rmassz(:)
-
! xmin
! if two chunks exclude this face for one of them
if(NCHUNKS == 1 .or. ichunk == CHUNK_AC) then
@@ -387,10 +450,18 @@
rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ if(SIMULATION_TYPE == 3 .and. ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + sngl(tx*weight)
+ b_rmassy(iglob) = b_rmassy(iglob) + sngl(ty*weight)
+ endif
else
rmassx(iglob) = rmassx(iglob) + tx*weight
rmassy(iglob) = rmassy(iglob) + ty*weight
rmassz(iglob) = rmassz(iglob) + tz*weight
+ if(SIMULATION_TYPE == 3 .and. ROTATION.and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + tx*weight
+ b_rmassy(iglob) = b_rmassy(iglob) + ty*weight
+ endif
endif
enddo
enddo
@@ -431,10 +502,18 @@
rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ if(SIMULATION_TYPE == 3 .and. ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + sngl(tx*weight)
+ b_rmassy(iglob) = b_rmassy(iglob) + sngl(ty*weight)
+ endif
else
rmassx(iglob) = rmassx(iglob) + tx*weight
rmassy(iglob) = rmassy(iglob) + ty*weight
rmassz(iglob) = rmassz(iglob) + tz*weight
+ if(SIMULATION_TYPE == 3 .and. ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + tx*weight
+ b_rmassy(iglob) = b_rmassy(iglob) + ty*weight
+ endif
endif
enddo
enddo
@@ -472,10 +551,18 @@
rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ if(SIMULATION_TYPE == 3 .and. ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + sngl(tx*weight)
+ b_rmassy(iglob) = b_rmassy(iglob) + sngl(ty*weight)
+ endif
else
rmassx(iglob) = rmassx(iglob) + tx*weight
rmassy(iglob) = rmassy(iglob) + ty*weight
rmassz(iglob) = rmassz(iglob) + tz*weight
+ if(SIMULATION_TYPE == 3 .and. ROTATION.and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + tx*weight
+ b_rmassy(iglob) = b_rmassy(iglob) + ty*weight
+ endif
endif
enddo
enddo
@@ -511,10 +598,18 @@
rmassx(iglob) = rmassx(iglob) + sngl(tx*weight)
rmassy(iglob) = rmassy(iglob) + sngl(ty*weight)
rmassz(iglob) = rmassz(iglob) + sngl(tz*weight)
+ if(SIMULATION_TYPE == 3 .and. ROTATION.and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + sngl(tx*weight)
+ b_rmassy(iglob) = b_rmassy(iglob) + sngl(ty*weight)
+ endif
else
rmassx(iglob) = rmassx(iglob) + tx*weight
rmassy(iglob) = rmassy(iglob) + ty*weight
rmassz(iglob) = rmassz(iglob) + tz*weight
+ if(SIMULATION_TYPE == 3 .and. ROTATION.and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ b_rmassx(iglob) = b_rmassx(iglob) + tx*weight
+ b_rmassy(iglob) = b_rmassy(iglob) + ty*weight
+ endif
endif
enddo
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -43,7 +43,7 @@
ner,ratio_sampling_array,doubling_index,r_bottom,r_top, &
this_region_has_a_doubling,ipass,ratio_divide_central_cube, &
CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,offset_proc_xi,offset_proc_eta,USE_FULL_TISO_MANTLE, &
- ATT1,ATT2,ATT3)
+ ATT1,ATT2,ATT3,SIMULATION_TYPE,USE_LDDRK,EXACT_MASS_MATRIX_FOR_ROTATION)
! creates the different regions of the mesh
@@ -163,6 +163,11 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
integer :: nglob_xy
+ ! mass matrices for backward simulation when SIMULATION_TYPE =3 and ROTATION is .true.
+ integer :: SIMULATION_TYPE,nglob_xy_backward
+ logical :: EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_rmassx,b_rmassy
+
! mass matrix and bathymetry for ocean load
integer nglob_oceans
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
@@ -977,8 +982,10 @@
!
! 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
+ nglob_xy = 1
+ nglob_xy_backward = 1
- if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+ if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS .and. (.not. USE_LDDRK)) then
select case(iregion_code)
case( IREGION_CRUST_MANTLE )
nglob_xy = nglob
@@ -989,6 +996,30 @@
nglob_xy = 1
endif
+ if(SIMULATION_TYPE /=3 .and. (.not. USE_LDDRK))then
+ if(ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ select case(iregion_code)
+ case( IREGION_CRUST_MANTLE,IREGION_INNER_CORE )
+ nglob_xy = nglob
+ case( IREGION_OUTER_CORE )
+ nglob_xy = 1
+ endselect
+ endif
+ endif
+
+ if(SIMULATION_TYPE ==3 .and. (.not. USE_LDDRK) )then
+ if(ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ select case(iregion_code)
+ case( IREGION_CRUST_MANTLE,IREGION_INNER_CORE )
+ nglob_xy = nglob
+ nglob_xy_backward = nglob
+ case( IREGION_OUTER_CORE )
+ nglob_xy = 1
+ nglob_xy_backward = 1
+ endselect
+ endif
+ endif
+
allocate(rmassx(nglob_xy),stat=ier)
if(ier /= 0) stop 'error in allocate 21'
allocate(rmassy(nglob_xy),stat=ier)
@@ -996,6 +1027,11 @@
allocate(rmassz(nglob),stat=ier)
if(ier /= 0) stop 'error in allocate 21'
+ allocate(b_rmassx(nglob_xy_backward),stat=ier)
+ if(ier /= 0) stop 'error in allocate b_21'
+ allocate(b_rmassy(nglob_xy_backward),stat=ier)
+ if(ier /= 0) stop 'error in allocate b_21'
+
! allocates ocean load mass matrix as well if oceans
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) then
nglob_oceans = nglob
@@ -1022,7 +1058,9 @@
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)
+ jacobian2D_bottom,jacobian2D_top,&
+ SIMULATION_TYPE,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ nglob_xy_backward,b_rmassx,b_rmassy)
! save the binary files
#ifdef USE_SERIAL_CASCADE_FOR_IOs
@@ -1051,7 +1089,9 @@
ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION, &
ATT1,ATT2,ATT3,size(tau_e_store,5),&
- NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso,myrank)
+ NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso,myrank,&
+ SIMULATION_TYPE,ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK,&
+ nglob_xy_backward,b_rmassx,b_rmassy)
#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)
@@ -1060,6 +1100,8 @@
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(b_rmassx,stat=ier); if(ier /= 0) stop 'error in deallocate b_rmassx'
+ deallocate(b_rmassy,stat=ier); if(ier /= 0) stop 'error in deallocate b_rmassy'
deallocate(rmass_ocean_load,stat=ier); if(ier /= 0) stop 'error in deallocate rmass_ocean_load'
! boundary mesh
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -712,7 +712,7 @@
this_region_has_a_doubling,ipass,ratio_divide_central_cube, &
CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA, &
mod(iproc_xi_slice(myrank),2),mod(iproc_eta_slice(myrank),2),USE_FULL_TISO_MANTLE, &
- ATT1,ATT2,ATT3)
+ ATT1,ATT2,ATT3,SIMULATION_TYPE,USE_LDDRK,EXACT_MASS_MATRIX_FOR_ROTATION)
enddo
! checks number of anisotropic elements found in the mantle
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -44,18 +44,21 @@
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE, &
ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION,ATT1,ATT2,ATT3,vnspec, &
- NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso,myrank)
+ NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso,myrank,&
+ SIMULATION_TYPE,ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK,&
+ nglob_xy_backward,b_rmassx,b_rmassy)
implicit none
include "constants.h"
- logical ATTENUATION
+ integer SIMULATION_TYPE
+ logical ATTENUATION,ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
character(len=150) prname
integer iregion_code,NCHUNKS
- integer nspec,nglob_xy,nglob,nspec_stacey,myrank
+ integer nspec,nglob_xy,nglob,nspec_stacey,myrank,nglob_xy_backward
integer npointot_oceans
! Stacey
@@ -97,6 +100,9 @@
real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
+! mass matrices for backward simulation when SIMULATION_TYPE =3 and ROTATION is .true.
+ real(kind=CUSTOM_REAL), dimension(nglob_xy_backward) :: b_rmassx,b_rmassy
+
! additional ocean load mass matrix
real(kind=CUSTOM_REAL) rmass_ocean_load(npointot_oceans)
@@ -257,13 +263,28 @@
!
! 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
+
+ if(.not. USE_LDDRK)then
+ if((NCHUNKS /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) .or. &
+ (ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. iregion_code == IREGION_CRUST_MANTLE) .or. &
+ (ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. iregion_code == IREGION_INNER_CORE)) then
+ write(27) rmassx
+ write(27) rmassy
+ endif
endif
write(27) rmassz
+ if(.not. USE_LDDRK)then
+ if(EXACT_MASS_MATRIX_FOR_ROTATION)then
+ if((SIMULATION_TYPE == 3 .and. (ROTATION .and. iregion_code == IREGION_CRUST_MANTLE)) .or. &
+ (SIMULATION_TYPE == 3 .and. (ROTATION .and. iregion_code == IREGION_INNER_CORE)))then
+ write(27) b_rmassx
+ write(27) b_rmassy
+ endif
+ endif
+ endif
+
! 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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -370,7 +370,7 @@
ibool_crust_mantle,ibelm_top_crust_mantle, &
updated_dof_ocean_load,NGLOB_XY, &
nspec_top, &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
implicit none
@@ -400,7 +400,7 @@
integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
logical, dimension(NGLOB_CRUST_MANTLE_OCEANS) :: updated_dof_ocean_load
- logical :: ABSORBING_CONDITIONS
+ logical :: ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
integer nspec_top
@@ -414,7 +414,8 @@
! initialize the updates
updated_dof_ocean_load(:) = .false.
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) .and. (.not. USE_LDDRK)) then
! for surface elements exactly at the top of the crust (ocean bottom)
do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -259,7 +259,7 @@
USE_LDDRK,INCREASE_CFL_FOR_LDDRK,ANISOTROPIC_KL,SAVE_TRANSVERSE_KL_ONLY,APPROXIMATE_HESS_KL, &
USE_FULL_TISO_MANTLE,SAVE_SOURCE_MASK,GPU_MODE,ADIOS_ENABLED,ADIOS_FOR_FORWARD_ARRAYS, &
ADIOS_FOR_MPI_ARRAYS,ADIOS_FOR_ARRAYS_SOLVER,ADIOS_FOR_AVS_DX,RATIO_BY_WHICH_TO_INCREASE_IT, &
- ATT1,ATT2,ATT3,ATT4,ATT5,RECOMPUTE_STRAIN_DO_NOT_STORE,EXACT_MASS_MATRIX_FOR_ROTATION)
+ ATT1,ATT2,ATT3,ATT4,ATT5,RECOMPUTE_STRAIN_DO_NOT_STORE,EXACT_MASS_MATRIX_FOR_ROTATION)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -796,9 +796,9 @@
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, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
endif
!-------------------------------------------------------------------------------------------------
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -7,7 +7,7 @@
! Newmark time scheme update
- do istage = 1, NSTAGE_TIME_SCHEME !ZN begin loop of istage
+ do istage = 1, NSTAGE_TIME_SCHEME ! begin loop of istage
if(USE_LDDRK)then
! mantle
accel_crust_mantle(:,:) = 0._CUSTOM_REAL
@@ -791,7 +791,8 @@
enddo
endif ! end of assembling forces with the central cube
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) .and. (.not. USE_LDDRK)) then
do i=1,NGLOB_CRUST_MANTLE
accel_crust_mantle(1,i) = accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
@@ -820,9 +821,9 @@
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, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
endif
!-------------------------------------------------------------------------------------------------
@@ -851,11 +852,19 @@
endif
! inner core
do i=1,NGLOB_INNER_CORE
- accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
- + two_omega_earth*veloc_inner_core(2,i)
- accel_inner_core(2,i) = accel_inner_core(2,i)*rmass_inner_core(i) &
- - two_omega_earth*veloc_inner_core(1,i)
- accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. .not. USE_LDDRK)then
+ accel_inner_core(1,i) = accel_inner_core(1,i)*rmassx_inner_core(i) &
+ + two_omega_earth*veloc_inner_core(2,i)
+ accel_inner_core(2,i) = accel_inner_core(2,i)*rmassy_inner_core(i) &
+ - two_omega_earth*veloc_inner_core(1,i)
+ accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ else
+ accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
+ + two_omega_earth*veloc_inner_core(2,i)
+ accel_inner_core(2,i) = accel_inner_core(2,i)*rmass_inner_core(i) &
+ - two_omega_earth*veloc_inner_core(1,i)
+ accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
+ endif
if(USE_LDDRK)then
veloc_inner_core_lddrk(:,i) = ALPHA_LDDRK(istage) * veloc_inner_core_lddrk(:,i) &
@@ -868,7 +877,7 @@
veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
endif
enddo
- enddo !ZN end loop of istage
+ enddo ! end loop of istage
! write the seismograms with time shift
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -792,7 +792,7 @@
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, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
ABSORBING_CONDITIONS)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -826,26 +826,44 @@
! ------------------- new non blocking implementation -------------------
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .and. .not. USE_LDDRK) then
- do i=1,NGLOB_CRUST_MANTLE
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassx_crust_mantle(i) &
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*b_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_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*b_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
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*b_rmassz_crust_mantle(i)
+ enddo
+ else
+ 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)*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
else
-
- do i=1,NGLOB_CRUST_MANTLE
- b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*rmassz_crust_mantle(i) &
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. .not. USE_LDDRK)then
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*b_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)*rmassz_crust_mantle(i) &
+ b_accel_crust_mantle(2,i) = b_accel_crust_mantle(2,i)*b_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
-
+ b_accel_crust_mantle(3,i) = b_accel_crust_mantle(3,i)*b_rmassz_crust_mantle(i)
+ enddo
+ else
+ do i=1,NGLOB_CRUST_MANTLE
+ b_accel_crust_mantle(1,i) = b_accel_crust_mantle(1,i)*b_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)*b_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)*b_rmassz_crust_mantle(i)
+ enddo
+ endif
endif
endif ! SIMULATION_TYPE == 3
@@ -854,12 +872,12 @@
if(OCEANS_VAL) then
if(SIMULATION_TYPE == 3) then
call compute_coupling_ocean(b_accel_crust_mantle, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,b_rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY, &
+ updated_dof_ocean_load,NGLOB_XY_CM, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- ABSORBING_CONDITIONS)
+ ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK)
endif
endif
@@ -879,11 +897,19 @@
enddo
! inner core
do i=1,NGLOB_INNER_CORE
- b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*rmass_inner_core(i) &
- + b_two_omega_earth*b_veloc_inner_core(2,i)
- b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*rmass_inner_core(i) &
- - b_two_omega_earth*b_veloc_inner_core(1,i)
- b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*rmass_inner_core(i)
+ if((ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) .and. .not. USE_LDDRK) then
+ b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*b_rmassx_inner_core(i) &
+ + b_two_omega_earth*b_veloc_inner_core(2,i)
+ b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*b_rmassy_inner_core(i) &
+ - b_two_omega_earth*b_veloc_inner_core(1,i)
+ b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*b_rmass_inner_core(i)
+ else
+ b_accel_inner_core(1,i) = b_accel_inner_core(1,i)*b_rmass_inner_core(i) &
+ + b_two_omega_earth*b_veloc_inner_core(2,i)
+ b_accel_inner_core(2,i) = b_accel_inner_core(2,i)*b_rmass_inner_core(i) &
+ - b_two_omega_earth*b_veloc_inner_core(1,i)
+ b_accel_inner_core(3,i) = b_accel_inner_core(3,i)*b_rmass_inner_core(i)
+ endif
b_veloc_inner_core(:,i) = b_veloc_inner_core(:,i) + b_deltatover2*b_accel_inner_core(:,i)
enddo
endif ! SIMULATION_TYPE == 3
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -44,26 +44,35 @@
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,NGLOB_XY,ABSORBING_CONDITIONS, &
- NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,npoin2D_max_all_CM_IC)
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS,NGLOB_XY_CM,ABSORBING_CONDITIONS, &
+ NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,npoin2D_max_all_CM_IC, &
+ SIMULATION_TYPE,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK,NGLOB_XY_CM_BACKWARD,&
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,&
+ NGLOB_XY_IC,rmassx_inner_core,rmassy_inner_core,&
+ NGLOB_XY_IC_BACKWARD,b_rmassx_inner_core,b_rmassy_inner_core)
use mpi
-
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
integer myrank,npoin2D_max_all_CM_IC
- integer NGLOB_XY
+ integer SIMULATION_TYPE,NGLOB_XY_CM,NGLOB_XY_CM_BACKWARD,NGLOB_XY_IC,NGLOB_XY_IC_BACKWARD
- logical ABSORBING_CONDITIONS
+ logical ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
- 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_XY_CM) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM_BACKWARD) :: b_rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM_BACKWARD) :: b_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_XY_IC) :: rmassx_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC) :: rmassy_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC_BACKWARD) :: b_rmassx_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC_BACKWARD) :: b_rmassy_inner_core
real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: rmass_inner_core
integer ichunk,iproc_xi,iproc_eta
@@ -140,7 +149,9 @@
endif
! crust and mantle
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) .and. NGLOB_CRUST_MANTLE > 0 &
+ .and. (.not. USE_LDDRK)) 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, &
@@ -184,6 +195,38 @@
NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE), &
NGLOB2DMAX_XY_CM_VAL,NCHUNKS_VAL)
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION &
+ .and. (.not. USE_LDDRK) .and. NGLOB_XY_CM_BACKWARD > 0)then
+ if(SIMULATION_TYPE == 3)then
+ call assemble_MPI_scalar_block(myrank,b_rmassx_crust_mantle,NGLOB_XY_CM_BACKWARD, &
+ 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_CM_VAL,NCHUNKS_VAL)
+ call assemble_MPI_scalar_block(myrank,b_rmassy_crust_mantle,NGLOB_XY_CM_BACKWARD, &
+ 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_CM_VAL,NCHUNKS_VAL)
+ endif
+ endif
+
! outer core
call assemble_MPI_scalar_block(myrank,rmass_outer_core,NGLOB_OUTER_CORE, &
iproc_xi,iproc_eta,ichunk,addressing, &
@@ -199,6 +242,70 @@
NGLOB2DMAX_XMIN_XMAX(IREGION_OUTER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_OUTER_CORE), &
NGLOB2DMAX_XY_OC_VAL,NCHUNKS_VAL)
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION &
+ .and. (.not. USE_LDDRK) .and. NGLOB_XY_IC > 0)then
+ ! inner core
+ call assemble_MPI_scalar_block(myrank,rmassx_inner_core,NGLOB_XY_IC, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi_inner_core,iboolright_xi_inner_core,iboolleft_eta_inner_core,iboolright_eta_inner_core, &
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
+ iboolfaces_inner_core,iboolcorner_inner_core, &
+ 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_INNER_CORE), &
+ NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE), &
+ NGLOB2DMAX_XY_IC_VAL,NCHUNKS_VAL)
+
+ call assemble_MPI_scalar_block(myrank,rmassy_inner_core,NGLOB_XY_IC, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi_inner_core,iboolright_xi_inner_core,iboolleft_eta_inner_core,iboolright_eta_inner_core, &
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
+ iboolfaces_inner_core,iboolcorner_inner_core, &
+ 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_INNER_CORE), &
+ NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE), &
+ NGLOB2DMAX_XY_IC_VAL,NCHUNKS_VAL)
+
+ if(SIMULATION_TYPE == 3 .and. NGLOB_XY_IC_BACKWARD > 0)then
+ call assemble_MPI_scalar_block(myrank,b_rmassx_inner_core,NGLOB_XY_IC_BACKWARD, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi_inner_core,iboolright_xi_inner_core,iboolleft_eta_inner_core,iboolright_eta_inner_core, &
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
+ iboolfaces_inner_core,iboolcorner_inner_core, &
+ 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_INNER_CORE), &
+ NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE), &
+ NGLOB2DMAX_XY_IC_VAL,NCHUNKS_VAL)
+
+ call assemble_MPI_scalar_block(myrank,b_rmassy_inner_core,NGLOB_XY_IC_BACKWARD, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi_inner_core,iboolright_xi_inner_core,iboolleft_eta_inner_core,iboolright_eta_inner_core, &
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
+ iboolfaces_inner_core,iboolcorner_inner_core, &
+ 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_INNER_CORE), &
+ NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE), &
+ NGLOB2DMAX_XY_IC_VAL,NCHUNKS_VAL)
+
+ endif
+
+ endif
+
! inner core
call assemble_MPI_scalar_block(myrank,rmass_inner_core,NGLOB_INNER_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 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -38,15 +38,17 @@
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)
+ ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS,LOCAL_PATH,ABSORBING_CONDITIONS,&
+ SIMULATION_TYPE,nglob_xy_backward,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ b_rmassx,b_rmassy)
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
- integer :: iregion_code,myrank
- integer :: nspec,nglob,nglob_xy
+ integer :: iregion_code,myrank,SIMULATION_TYPE
+ integer :: nspec,nglob,nglob_xy,nglob_xy_backward
integer :: nspec_iso,nspec_tiso,nspec_ani
! Stacey
@@ -78,9 +80,12 @@
! mass matrices and additional ocean load mass matrix
real(kind=CUSTOM_REAL), dimension(nglob_xy) :: rmassx,rmassy
+ real(kind=CUSTOM_REAL), dimension(nglob_xy_backward) :: b_rmassx,b_rmassy
real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
real(kind=CUSTOM_REAL), dimension(nglob) :: rmass_ocean_load
+ logical :: EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
+
! flags to know if we should read Vs and anisotropy arrays
logical :: READ_KAPPA_MU,READ_TISO,ABSORBING_CONDITIONS,TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS
@@ -175,13 +180,32 @@
!
! 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
+! if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) then
+! read(IIN) rmassx
+! read(IIN) rmassy
+! endif
+
+ if(.not. USE_LDDRK)then
+ if((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. iregion_code == IREGION_CRUST_MANTLE) .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. iregion_code == IREGION_INNER_CORE)) then
+ read(IIN) rmassx
+ read(IIN) rmassy
+ endif
endif
read(IIN) rmassz
+ if(.not. USE_LDDRK)then
+ if((SIMULATION_TYPE == 3 .and. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. iregion_code == IREGION_CRUST_MANTLE)) .or. &
+ (SIMULATION_TYPE == 3 .and. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. iregion_code == IREGION_INNER_CORE)))then
+ read(IIN) b_rmassx
+ read(IIN) b_rmassy
+ endif
+ endif
+
! 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 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -59,14 +59,18 @@
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,NGLOB_XY)
+ ABSORBING_CONDITIONS,LOCAL_PATH,NGLOB_XY_CM,&
+ SIMULATION_TYPE,NGLOB_XY_CM_BACKWARD,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,&
+ NGLOB_XY_IC,rmassx_inner_core,rmassy_inner_core,&
+ NGLOB_XY_IC_BACKWARD,b_rmassx_inner_core,b_rmassy_inner_core)
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
- integer myrank
+ integer myrank,SIMULATION_TYPE
! Stacey
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STACEY) :: &
@@ -108,10 +112,12 @@
!
! 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
+ integer :: NGLOB_DUMMY,NGLOB_XY_CM,NGLOB_XY_CM_BACKWARD
+ logical :: EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM) :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM_BACKWARD) :: b_rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_CM_BACKWARD) :: b_rmassy_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
! additional mass matrix for ocean load
@@ -148,6 +154,11 @@
integer, dimension(NSPEC_INNER_CORE) :: idoubling_inner_core
logical, dimension(NSPEC_INNER_CORE) :: ispec_is_tiso_inner_core
+ integer :: NGLOB_XY_IC,NGLOB_XY_IC_BACKWARD
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC) :: rmassx_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC) :: rmassy_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC_BACKWARD) :: b_rmassx_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_XY_IC_BACKWARD) :: b_rmassy_inner_core
real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: rmass_inner_core
logical ABSORBING_CONDITIONS
@@ -209,11 +220,14 @@
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,ispec_is_tiso_crust_mantle,NGLOB_XY,NGLOB_CRUST_MANTLE, &
+ ibool_crust_mantle,dummy_i,ispec_is_tiso_crust_mantle,NGLOB_XY_CM,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)
+ ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS,&
+ SIMULATION_TYPE,NGLOB_XY_CM_BACKWARD,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle)
+
! synchronizes processes
call sync_all()
@@ -244,10 +258,12 @@
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,NGLOB_XY,NGLOB_OUTER_CORE, &
+ ibool_outer_core,idoubling_outer_core,ispec_is_tiso_outer_core,1,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)
+ ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS,&
+ SIMULATION_TYPE,NGLOB_DUMMY,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ dummy_rmass,dummy_rmass)
! synchronizes processes
call sync_all()
@@ -280,10 +296,12 @@
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,NGLOB_XY,NGLOB_INNER_CORE, &
- dummy_rmass,dummy_rmass,rmass_inner_core,rmass_ocean_load,NSPEC_INNER_CORE, &
+ ibool_inner_core,idoubling_inner_core,ispec_is_tiso_inner_core,NGLOB_XY_IC,NGLOB_INNER_CORE, &
+ rmassx_inner_core,rmassy_inner_core,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)
+ ANISOTROPIC_INNER_CORE_VAL,OCEANS_VAL,LOCAL_PATH,ABSORBING_CONDITIONS,&
+ SIMULATION_TYPE,NGLOB_XY_IC_BACKWARD,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ b_rmassx_inner_core,b_rmassy_inner_core)
! check that the number of points in this slice is correct
if(minval(ibool_crust_mantle(:,:,:,:)) /= 1 .or. &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2013-07-12 11:12:06 UTC (rev 22573)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2013-07-12 17:09:21 UTC (rev 22574)
@@ -555,11 +555,21 @@
!
! 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
+!
+! in the case of ROTATION, we should add two_omega_earth*deltat/2 contribution to rmassx & rmassy
+! thus in this case rmassx & rmassy will be used
+!
+! in the case of ROTAION and SIMULATION_TYPE == 3, we should add b_two_omega_earth*deltat/2 contribution to &
+! b_rmassx & b_rmassy
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx_crust_mantle
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassy_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_rmassy_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: b_rmassz_crust_mantle
+ equivalence(rmassz_crust_mantle,b_rmassz_crust_mantle)
- integer :: NGLOB_XY
+ integer :: NGLOB_XY_CM,NGLOB_XY_CM_BACKWARD
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
@@ -614,7 +624,14 @@
logical, dimension(NSPEC_INNER_CORE) :: ispec_is_tiso_inner_core ! only needed for computer_boundary_kernel() routine
! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: rmass_inner_core
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx_inner_core
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassy_inner_core
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_rmassx_inner_core
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_rmassy_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: rmass_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: b_rmass_inner_core
+ integer :: NGLOB_XY_IC,NGLOB_XY_IC_BACKWARD
+ equivalence(rmass_inner_core,b_rmass_inner_core)
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: &
@@ -1179,17 +1196,54 @@
!
! 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
+
+ NGLOB_XY_CM = 1
+ NGLOB_XY_IC = 1
+ NGLOB_XY_CM_BACKWARD = 1
+ NGLOB_XY_IC_BACKWARD = 1
+
+ if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS .and. (.not. USE_LDDRK)) then
+ NGLOB_XY_CM = NGLOB_CRUST_MANTLE
else
- NGLOB_XY = 1
+ NGLOB_XY_CM = 1
endif
- allocate(rmassx_crust_mantle(NGLOB_XY),stat=ier)
+ if(SIMULATION_TYPE /=3 .and. (.not. USE_LDDRK) .and. EXACT_MASS_MATRIX_FOR_ROTATION)then
+ if(ROTATION_VAL)then
+ NGLOB_XY_CM = NGLOB_CRUST_MANTLE
+ NGLOB_XY_IC = NGLOB_INNER_CORE
+ endif
+ endif
+
+ if(SIMULATION_TYPE ==3 .and. (.not. USE_LDDRK) .and. EXACT_MASS_MATRIX_FOR_ROTATION )then
+ if(ROTATION_VAL)then
+ NGLOB_XY_CM = NGLOB_CRUST_MANTLE
+ NGLOB_XY_IC = NGLOB_INNER_CORE
+ NGLOB_XY_CM_BACKWARD = NGLOB_CRUST_MANTLE
+ NGLOB_XY_IC_BACKWARD = NGLOB_INNER_CORE
+ endif
+ endif
+
+ allocate(rmassx_crust_mantle(NGLOB_XY_CM),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating rmassx_crust_mantle')
- allocate(rmassy_crust_mantle(NGLOB_XY),stat=ier)
+ allocate(rmassy_crust_mantle(NGLOB_XY_CM),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating rmassy_crust_mantle')
+ allocate(b_rmassx_crust_mantle(NGLOB_XY_CM_BACKWARD),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_rmassx_crust_mantle')
+ allocate(b_rmassy_crust_mantle(NGLOB_XY_CM_BACKWARD),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_rmassy_crust_mantle')
+
+ allocate(rmassx_inner_core(NGLOB_XY_IC),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating rmassx_inner_core')
+ allocate(rmassy_inner_core(NGLOB_XY_IC),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating rmassy_inner_core')
+
+ allocate(b_rmassx_inner_core(NGLOB_XY_IC_BACKWARD),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_rmassx_inner_core')
+ allocate(b_rmassy_inner_core(NGLOB_XY_IC_BACKWARD),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_rmassy_inner_core')
+
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, &
@@ -1224,7 +1278,11 @@
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,NGLOB_XY)
+ ABSORBING_CONDITIONS,LOCAL_PATH,NGLOB_XY_CM,&
+ SIMULATION_TYPE,NGLOB_XY_CM_BACKWARD,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,&
+ NGLOB_XY_IC,rmassx_inner_core,rmassy_inner_core,&
+ NGLOB_XY_IC_BACKWARD,b_rmassx_inner_core,b_rmassy_inner_core)
! read 2-D addressing for summation between slices with MPI
call read_mesh_databases_addressing(myrank, &
@@ -1744,8 +1802,12 @@
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,NGLOB_XY,ABSORBING_CONDITIONS, &
- NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,npoin2D_max_all_CM_IC)
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS,NGLOB_XY_CM,ABSORBING_CONDITIONS, &
+ NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,npoin2D_max_all_CM_IC,&
+ SIMULATION_TYPE,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK,NGLOB_XY_CM_BACKWARD,&
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,&
+ NGLOB_XY_IC,rmassx_inner_core,rmassy_inner_core,&
+ NGLOB_XY_IC_BACKWARD,b_rmassx_inner_core,b_rmassy_inner_core)
! mass matrix including central cube
if(INCLUDE_CENTRAL_CUBE) then
@@ -1797,6 +1859,69 @@
sender_from_slices_to_cube,ibool_central_cube, &
buffer_slices,buffer_slices2,buffer_all_cube_from_slices)
+ if(ROTATION_VAL .and. (.not. USE_LDDRK) .and. EXACT_MASS_MATRIX_FOR_ROTATION &
+ .and. NGLOB_XY_IC > 0)then
+
+ call prepare_timerun_centralcube(myrank,rmassx_inner_core, &
+ iproc_xi,iproc_eta,ichunk, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM, &
+ addressing,ibool_inner_core,idoubling_inner_core, &
+ xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+ nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
+ nspec2D_ymin_inner_core,nspec2D_ymax_inner_core, &
+ ibelm_xmin_inner_core,ibelm_xmax_inner_core, &
+ ibelm_ymin_inner_core,ibelm_ymax_inner_core,ibelm_bottom_inner_core, &
+ nb_msgs_theor_in_cube,non_zero_nb_msgs_theor_in_cube, &
+ npoin2D_cube_from_slices,receiver_cube_from_slices, &
+ sender_from_slices_to_cube,ibool_central_cube, &
+ buffer_slices,buffer_slices2,buffer_all_cube_from_slices)
+
+ call prepare_timerun_centralcube(myrank,rmassy_inner_core, &
+ iproc_xi,iproc_eta,ichunk, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM, &
+ addressing,ibool_inner_core,idoubling_inner_core, &
+ xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+ nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
+ nspec2D_ymin_inner_core,nspec2D_ymax_inner_core, &
+ ibelm_xmin_inner_core,ibelm_xmax_inner_core, &
+ ibelm_ymin_inner_core,ibelm_ymax_inner_core,ibelm_bottom_inner_core, &
+ nb_msgs_theor_in_cube,non_zero_nb_msgs_theor_in_cube, &
+ npoin2D_cube_from_slices,receiver_cube_from_slices, &
+ sender_from_slices_to_cube,ibool_central_cube, &
+ buffer_slices,buffer_slices2,buffer_all_cube_from_slices)
+
+ if(SIMULATION_TYPE == 3 .and. NGLOB_XY_IC_BACKWARD > 0)then
+
+ call prepare_timerun_centralcube(myrank,b_rmassx_inner_core, &
+ iproc_xi,iproc_eta,ichunk, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM, &
+ addressing,ibool_inner_core,idoubling_inner_core, &
+ xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+ nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
+ nspec2D_ymin_inner_core,nspec2D_ymax_inner_core, &
+ ibelm_xmin_inner_core,ibelm_xmax_inner_core, &
+ ibelm_ymin_inner_core,ibelm_ymax_inner_core,ibelm_bottom_inner_core, &
+ nb_msgs_theor_in_cube,non_zero_nb_msgs_theor_in_cube, &
+ npoin2D_cube_from_slices,receiver_cube_from_slices, &
+ sender_from_slices_to_cube,ibool_central_cube, &
+ buffer_slices,buffer_slices2,buffer_all_cube_from_slices)
+
+ call prepare_timerun_centralcube(myrank,b_rmassy_inner_core, &
+ iproc_xi,iproc_eta,ichunk, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM, &
+ addressing,ibool_inner_core,idoubling_inner_core, &
+ xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+ nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
+ nspec2D_ymin_inner_core,nspec2D_ymax_inner_core, &
+ ibelm_xmin_inner_core,ibelm_xmax_inner_core, &
+ ibelm_ymin_inner_core,ibelm_ymax_inner_core,ibelm_bottom_inner_core, &
+ nb_msgs_theor_in_cube,non_zero_nb_msgs_theor_in_cube, &
+ npoin2D_cube_from_slices,receiver_cube_from_slices, &
+ sender_from_slices_to_cube,ibool_central_cube, &
+ buffer_slices,buffer_slices2,buffer_all_cube_from_slices)
+ endif
+ endif
+
call fix_non_blocking_central_cube(is_on_a_slice_edge_inner_core, &
ibool_inner_core,NSPEC_INNER_CORE,NGLOB_INNER_CORE,nb_msgs_theor_in_cube,ibelm_bottom_inner_core, &
idoubling_inner_core,npoin2D_cube_from_slices,ibool_central_cube, &
@@ -1824,7 +1949,9 @@
endif
! add C*delta/2 contribution to the mass matrices on the Stacey edges
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if(((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) &
+ .and. (.not. USE_LDDRK)) 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) &
@@ -1842,11 +1969,42 @@
if(OCEANS_VAL) rmass_ocean_load = 1._CUSTOM_REAL / rmass_ocean_load
! add C*delta/2 contribution to the mass matrices on the Stacey edges
- if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
+ if(((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) &
+ .and. (.not. USE_LDDRK)) then
rmassx_crust_mantle = 1._CUSTOM_REAL / rmassx_crust_mantle
rmassy_crust_mantle = 1._CUSTOM_REAL / rmassy_crust_mantle
endif
+ if(.not. USE_LDDRK)then
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION .and. NGLOB_XY_IC > 0) then
+ if(minval(rmassx_inner_core) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the rmassx_inner_core')
+ if(minval(rmassy_inner_core) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the rmassy_inner_core')
+ rmassx_inner_core = 1._CUSTOM_REAL / rmassx_inner_core
+ rmassy_inner_core = 1._CUSTOM_REAL / rmassy_inner_core
+ if(SIMULATION_TYPE == 3 .and. NGLOB_XY_IC_BACKWARD > 0)then
+ if(minval(b_rmassx_inner_core) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the b_rmassx_inner_core')
+ if(minval(b_rmassy_inner_core) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the b_rmassy_inner_core')
+ b_rmassx_inner_core = 1._CUSTOM_REAL / b_rmassx_inner_core
+ b_rmassy_inner_core = 1._CUSTOM_REAL / b_rmassy_inner_core
+ endif
+ endif
+
+ if(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION &
+ .and. (.not. USE_LDDRK) .and. SIMULATION_TYPE == 3 .and. NGLOB_XY_CM_BACKWARD > 0)then
+ if(minval(b_rmassx_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the b_crust_mantle')
+ if(minval(b_rmassy_crust_mantle) <= 0._CUSTOM_REAL) &
+ call exit_MPI(myrank,'negative mass matrix term for the b_crust_mantle')
+ b_rmassx_crust_mantle = 1._CUSTOM_REAL / b_rmassx_crust_mantle
+ b_rmassy_crust_mantle = 1._CUSTOM_REAL / b_rmassy_crust_mantle
+ endif
+ 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
@@ -2332,7 +2490,7 @@
!! DK DK this first part handles the cases SIMULATION_TYPE == 1 and SIMULATION_TYPE == 2
!! DK DK it also handles the cases NOISE_TOMOGRAPHY == 1 and NOISE_TOMOGRAPHY == 2
!! DK DK
- if(USE_LDDRK)then
+ if(USE_LDDRK .or. EXACT_MASS_MATRIX_FOR_ROTATION)then
include "part1_undo_att.f90"
else
include "part1_classical.f90"
More information about the CIG-COMMITS
mailing list