[cig-commits] r15043 - seismo/3D/SPECFEM3D_SESAME/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon May 25 09:12:49 PDT 2009
Author: dkomati1
Date: 2009-05-25 09:12:48 -0700 (Mon, 25 May 2009)
New Revision: 15043
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
Removed:
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_scalar.f90
seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_vector.f90
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90
seismo/3D/SPECFEM3D_SESAME/trunk/configure
seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac
seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
seismo/3D/SPECFEM3D_SESAME/trunk/create_header_file.f90
seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
Log:
removed old routines related to a regular mesh;
added an option to turn the optimized products of Deville on or off.
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-05-25 16:12:48 UTC (rev 15043)
@@ -109,7 +109,8 @@
$O/write_AVS_DX_surface_data.o \
$O/write_seismograms.o \
$O/compute_boundary_kernel.o \
- $O/compute_forces.o \
+ $O/compute_forces_no_Deville.o \
+ $O/compute_forces_with_Deville.o \
$(EMPTY_MACRO)
# solver objects with statically allocated arrays; dependent upon
@@ -427,9 +428,12 @@
$O/compute_boundary_kernel.o: constants.h compute_boundary_kernel.f90
${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o compute_boundary_kernel.f90
-$O/compute_forces.o: constants.h compute_forces.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces.o compute_forces.f90
+$O/compute_forces_no_Deville.o: constants.h compute_forces_no_Deville.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_no_Deville.o compute_forces_no_Deville.f90
+$O/compute_forces_with_Deville.o: constants.h compute_forces_with_Deville.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_with_Deville.o compute_forces_with_Deville.f90
+
$O/combine_vol_data.o: constants.h combine_vol_data.f90
${FCCOMPILE_CHECK} -c -o $O/combine_vol_data.o combine_vol_data.f90
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_scalar.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_scalar.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_scalar.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -24,188 +24,14 @@
!=====================================================================
!----
-!---- assemble the contributions between slices and chunks using MPI
+!---- assemble the contributions between slices using non-blocking MPI
!----
- subroutine assemble_MPI_scalar(array_val,iproc_xi,iproc_eta,addressing, &
- iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
- buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_xi,npoin2D_eta, &
- NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
-
- implicit none
-
- include "constants.h"
-
-! include values created by the mesher
- include "OUTPUT_FILES/values_from_mesher.h"
-
-! array to assemble
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB_VAL) :: array_val
-
- integer iproc_xi,iproc_eta
- integer npoin2D_xi,npoin2D_eta
-
- integer NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY
- integer NPROC_XI,NPROC_ETA
-
-! for addressing of the slices
- integer, dimension(0:NPROC_XI-1,0:NPROC_ETA-1) :: addressing
-
-! 2-D addressing and buffers for summation between slices
- integer, dimension(NPOIN2DMAX_XMIN_XMAX) :: iboolleft_xi,iboolright_xi
- integer, dimension(NPOIN2DMAX_YMIN_YMAX) :: iboolleft_eta,iboolright_eta
-
- real(kind=CUSTOM_REAL), dimension(NPOIN2DMAX_XY) :: buffer_send_faces_scalar,buffer_received_faces_scalar
-
- integer ipoin
- integer sender,receiver
-
- integer, external :: proc_null
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-! here we have to assemble all the contributions between slices using MPI
-
-!----
-!---- first assemble along xi using the 2-D topology
-!----
-
-! assemble along xi only if more than one slice
- if(NPROC_XI > 1) then
-
-! slices copy the right face into the buffer
- do ipoin=1,npoin2D_xi
- buffer_send_faces_scalar(ipoin) = array_val(iboolright_xi(ipoin))
- enddo
-
-! send messages forward along each row
- if(iproc_xi == 0) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi - 1,iproc_eta)
- endif
- if(iproc_xi == NPROC_XI-1) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi + 1,iproc_eta)
- endif
- call sendrecv_all_cr(buffer_send_faces_scalar,npoin2D_xi,receiver,itag2, &
- buffer_received_faces_scalar,npoin2D_xi,sender,itag)
-
-! all slices add the buffer received to the contributions on the left face
- if(iproc_xi > 0) then
- do ipoin=1,npoin2D_xi
- array_val(iboolleft_xi(ipoin)) = array_val(iboolleft_xi(ipoin)) + &
- buffer_received_faces_scalar(ipoin)
- enddo
- endif
-
-! the contributions are correctly assembled on the left side of each slice
-! now we have to send the result back to the sender
-! all slices copy the left face into the buffer
- do ipoin=1,npoin2D_xi
- buffer_send_faces_scalar(ipoin) = array_val(iboolleft_xi(ipoin))
- enddo
-
-! send messages backward along each row
- if(iproc_xi == NPROC_XI-1) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi + 1,iproc_eta)
- endif
- if(iproc_xi == 0) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi - 1,iproc_eta)
- endif
- call sendrecv_all_cr(buffer_send_faces_scalar,npoin2D_xi,receiver,itag2, &
- buffer_received_faces_scalar,npoin2D_xi,sender,itag)
-
-! all slices copy the buffer received to the contributions on the right face
- if(iproc_xi < NPROC_XI-1) then
- do ipoin=1,npoin2D_xi
- array_val(iboolright_xi(ipoin)) = buffer_received_faces_scalar(ipoin)
- enddo
- endif
-
- endif
-
-!----
-!---- then assemble along eta using the 2-D topology
-!----
-
-! assemble along eta only if more than one slice
- if(NPROC_ETA > 1) then
-
-! slices copy the right face into the buffer
- do ipoin=1,npoin2D_eta
- buffer_send_faces_scalar(ipoin) = array_val(iboolright_eta(ipoin))
- enddo
-
-! send messages forward along each row
- if(iproc_eta == 0) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi,iproc_eta - 1)
- endif
- if(iproc_eta == NPROC_ETA-1) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi,iproc_eta + 1)
- endif
- call sendrecv_all_cr(buffer_send_faces_scalar,npoin2D_eta,receiver,itag2, &
- buffer_received_faces_scalar,npoin2D_eta,sender,itag)
-
-! all slices add the buffer received to the contributions on the left face
- if(iproc_eta > 0) then
- do ipoin=1,npoin2D_eta
- array_val(iboolleft_eta(ipoin)) = array_val(iboolleft_eta(ipoin)) + &
- buffer_received_faces_scalar(ipoin)
- enddo
- endif
-
-! the contributions are correctly assembled on the left side of each slice
-! now we have to send the result back to the sender
-! all slices copy the left face into the buffer
- do ipoin=1,npoin2D_eta
- buffer_send_faces_scalar(ipoin) = array_val(iboolleft_eta(ipoin))
- enddo
-
-! send messages backward along each row
- if(iproc_eta == NPROC_ETA-1) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi,iproc_eta + 1)
- endif
- if(iproc_eta == 0) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi,iproc_eta - 1)
- endif
- call sendrecv_all_cr(buffer_send_faces_scalar,npoin2D_eta,receiver,itag2, &
- buffer_received_faces_scalar,npoin2D_eta,sender,itag)
-
-! all slices copy the buffer received to the contributions on the right face
- if(iproc_eta < NPROC_ETA-1) then
- do ipoin=1,npoin2D_eta
- array_val(iboolright_eta(ipoin)) = buffer_received_faces_scalar(ipoin)
- enddo
- endif
-
- endif
-
- end subroutine assemble_MPI_scalar
-
-!
-!----
-!
-
subroutine assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,array_val, &
buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh, &
ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
- request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
- )
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
implicit none
@@ -359,7 +185,7 @@
call wait_req(request_send_scalar_ext_mesh(iinterface))
enddo
-
endif
end subroutine assemble_MPI_scalar_i_ext_mesh
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_vector.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_vector.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/assemble_MPI_vector.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -24,188 +24,14 @@
!=====================================================================
!----
-!---- assemble the contributions between slices and chunks using MPI
+!---- assemble the contributions between slices using non-blocking MPI
!----
- subroutine assemble_MPI_vector(array_val,iproc_xi,iproc_eta,addressing, &
- iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
- buffer_send_faces_vector,buffer_received_faces_vector,npoin2D_xi,npoin2D_eta, &
- NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
-
- implicit none
-
- include "constants.h"
-
-! include values created by the mesher
- include "OUTPUT_FILES/values_from_mesher.h"
-
-! array to assemble
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB_VAL) :: array_val
-
- integer iproc_xi,iproc_eta
- integer npoin2D_xi,npoin2D_eta
-
- integer NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY
- integer NPROC_XI,NPROC_ETA
-
-! for addressing of the slices
- integer, dimension(0:NPROC_XI-1,0:NPROC_ETA-1) :: addressing
-
-! 2-D addressing and buffers for summation between slices
- integer, dimension(NPOIN2DMAX_XMIN_XMAX) :: iboolleft_xi,iboolright_xi
- integer, dimension(NPOIN2DMAX_YMIN_YMAX) :: iboolleft_eta,iboolright_eta
-
- real(kind=CUSTOM_REAL), dimension(NDIM,NPOIN2DMAX_XY) :: buffer_send_faces_vector,buffer_received_faces_vector
-
- integer ipoin
- integer sender,receiver
-
- integer, external :: proc_null
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-! here we have to assemble all the contributions between slices using MPI
-
-!----
-!---- first assemble along xi using the 2-D topology
-!----
-
-! assemble along xi only if more than one slice
- if(NPROC_XI > 1) then
-
-! slices copy the right face into the buffer
- do ipoin=1,npoin2D_xi
- buffer_send_faces_vector(:,ipoin) = array_val(:,iboolright_xi(ipoin))
- enddo
-
-! send messages forward along each row
- if(iproc_xi == 0) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi - 1,iproc_eta)
- endif
- if(iproc_xi == NPROC_XI-1) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi + 1,iproc_eta)
- endif
- call sendrecv_all_cr(buffer_send_faces_vector,NDIM*npoin2D_xi,receiver,itag2, &
- buffer_received_faces_vector,NDIM*npoin2D_xi,sender,itag)
-
-! all slices add the buffer received to the contributions on the left face
- if(iproc_xi > 0) then
- do ipoin=1,npoin2D_xi
- array_val(:,iboolleft_xi(ipoin)) = array_val(:,iboolleft_xi(ipoin)) + &
- buffer_received_faces_vector(:,ipoin)
- enddo
- endif
-
-! the contributions are correctly assembled on the left side of each slice
-! now we have to send the result back to the sender
-! all slices copy the left face into the buffer
- do ipoin=1,npoin2D_xi
- buffer_send_faces_vector(:,ipoin) = array_val(:,iboolleft_xi(ipoin))
- enddo
-
-! send messages backward along each row
- if(iproc_xi == NPROC_XI-1) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi + 1,iproc_eta)
- endif
- if(iproc_xi == 0) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi - 1,iproc_eta)
- endif
- call sendrecv_all_cr(buffer_send_faces_vector,NDIM*npoin2D_xi,receiver,itag2, &
- buffer_received_faces_vector,NDIM*npoin2D_xi,sender,itag)
-
-! all slices copy the buffer received to the contributions on the right face
- if(iproc_xi < NPROC_XI-1) then
- do ipoin=1,npoin2D_xi
- array_val(:,iboolright_xi(ipoin)) = buffer_received_faces_vector(:,ipoin)
- enddo
- endif
-
- endif
-
-!----
-!---- then assemble along eta using the 2-D topology
-!----
-
-! assemble along eta only if more than one slice
- if(NPROC_ETA > 1) then
-
-! slices copy the right face into the buffer
- do ipoin=1,npoin2D_eta
- buffer_send_faces_vector(:,ipoin) = array_val(:,iboolright_eta(ipoin))
- enddo
-
-! send messages forward along each row
- if(iproc_eta == 0) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi,iproc_eta - 1)
- endif
- if(iproc_eta == NPROC_ETA-1) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi,iproc_eta + 1)
- endif
- call sendrecv_all_cr(buffer_send_faces_vector,NDIM*npoin2D_eta,receiver,itag2, &
- buffer_received_faces_vector,NDIM*npoin2D_eta,sender,itag)
-
-! all slices add the buffer received to the contributions on the left face
- if(iproc_eta > 0) then
- do ipoin=1,npoin2D_eta
- array_val(:,iboolleft_eta(ipoin)) = array_val(:,iboolleft_eta(ipoin)) + &
- buffer_received_faces_vector(:,ipoin)
- enddo
- endif
-
-! the contributions are correctly assembled on the left side of each slice
-! now we have to send the result back to the sender
-! all slices copy the left face into the buffer
- do ipoin=1,npoin2D_eta
- buffer_send_faces_vector(:,ipoin) = array_val(:,iboolleft_eta(ipoin))
- enddo
-
-! send messages backward along each row
- if(iproc_eta == NPROC_ETA-1) then
- sender = proc_null()
- else
- sender = addressing(iproc_xi,iproc_eta + 1)
- endif
- if(iproc_eta == 0) then
- receiver = proc_null()
- else
- receiver = addressing(iproc_xi,iproc_eta - 1)
- endif
- call sendrecv_all_cr(buffer_send_faces_vector,NDIM*npoin2D_eta,receiver,itag2, &
- buffer_received_faces_vector,NDIM*npoin2D_eta,sender,itag)
-
-! all slices copy the buffer received to the contributions on the right face
- if(iproc_eta < NPROC_ETA-1) then
- do ipoin=1,npoin2D_eta
- array_val(:,iboolright_eta(ipoin)) = buffer_received_faces_vector(:,ipoin)
- enddo
- endif
-
- endif
-
- end subroutine assemble_MPI_vector
-
-!
-!----
-!
-
subroutine assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,array_val, &
buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
- request_send_vector_ext_mesh,request_recv_vector_ext_mesh &
- )
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
implicit none
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -1,383 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-subroutine compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- kappastore,mustore,jacobian,ibool,ispec_is_inner,phase_is_inner, &
- NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
-
- implicit none
-
- include "constants.h"
-
- integer :: NSPEC_AB,NGLOB_AB
-
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
-
-! arrays with mesh parameters per slice
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
- kappastore,mustore,jacobian
-
-! array with derivatives of Lagrange polynomials and precalculated products
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT
- 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), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
- newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
-
-! communication overlap
- logical, dimension(NSPEC_AB) :: ispec_is_inner
- logical :: phase_is_inner
-
-! source
- integer :: NSOURCES,myrank,it
- integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
- double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
- double precision, dimension(3,3,NSOURCES) :: nu_source
- double precision, dimension(NSOURCES) :: hdur
- double precision :: dt
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
-
- integer ispec,iglob
- integer i,j,k
-
- real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
- real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
-
- real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
- real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
-
- real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
-
- real(kind=CUSTOM_REAL) fac1,fac2,fac3
-
- real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
- real(kind=CUSTOM_REAL) kappal
-
- integer :: isource
- double precision :: t0,f0
-
-
- do ispec = 1,NSPEC_AB
-
- if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
- dummyx_loc(i,j,k) = displ(1,iglob)
- dummyy_loc(i,j,k) = displ(2,iglob)
- dummyz_loc(i,j,k) = displ(3,iglob)
- enddo
- enddo
- enddo
-
-!! DK DK subroutines adapted from Deville, Fischer and Mund, High-order methods
-!! DK DK for incompressible fluid flow, Cambridge University Press (2002),
-!! DK DK pages 386 and 389 and Figure 8.3.1
- call mxm_m1_m2_5points(hprime_xx,dummyx_loc,dummyy_loc,dummyz_loc,tempx1,tempy1,tempz1)
-
- do k = 1,NGLLX
- call mxm_m1_m1_5points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
- hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
- enddo
-
- call mxm_m2_m1_5points(dummyx_loc,dummyy_loc,dummyz_loc,hprime_xxT,tempx3,tempy3,tempz3)
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
-! get derivatives of ux, uy and uz with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
-
- duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
- duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
- duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
-
- duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
- duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
- duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
-
- duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
- duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
- duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
-
-! precompute some sums to save CPU time
- duxdxl_plus_duydyl = duxdxl + duydyl
- duxdxl_plus_duzdzl = duxdxl + duzdzl
- duydyl_plus_duzdzl = duydyl + duzdzl
- duxdyl_plus_duydxl = duxdyl + duydxl
- duzdxl_plus_duxdzl = duzdxl + duxdzl
- duzdyl_plus_duydzl = duzdyl + duydzl
-
- kappal = kappastore(i,j,k,ispec)
- mul = mustore(i,j,k,ispec)
-
- lambdalplus2mul = kappal + FOUR_THIRDS * mul
- lambdal = lambdalplus2mul - 2.*mul
-
-! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
- sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
- sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
-
- sigma_xy = mul*duxdyl_plus_duydxl
- sigma_xz = mul*duzdxl_plus_duxdzl
- sigma_yz = mul*duzdyl_plus_duydzl
-
-! form dot product with test vector, symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
-
- enddo
- enddo
- enddo
-
-!! DK DK subroutines adapted from Deville, Fischer and Mund, High-order methods
-!! DK DK for incompressible fluid flow, Cambridge University Press (2002),
-!! DK DK pages 386 and 389 and Figure 8.3.1
- call mxm_m1_m2_5points(hprimewgll_xxT,tempx1,tempy1,tempz1,newtempx1,newtempy1,newtempz1)
-
- do k = 1,NGLLX
- call mxm_m1_m1_5points(tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k), &
- hprimewgll_xx,newtempx2(1,1,k),newtempy2(1,1,k),newtempz2(1,1,k))
- enddo
-
- call mxm_m2_m1_5points(tempx3,tempy3,tempz3,hprimewgll_xx,newtempx3,newtempy3,newtempz3)
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
-
-! sum contributions from each element to the global mesh using indirect addressing
- iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = accel(1,iglob) - fac1*newtempx1(i,j,k) - fac2*newtempx2(i,j,k) - fac3*newtempx3(i,j,k)
- accel(2,iglob) = accel(2,iglob) - fac1*newtempy1(i,j,k) - fac2*newtempy2(i,j,k) - fac3*newtempy3(i,j,k)
- accel(3,iglob) = accel(3,iglob) - fac1*newtempz1(i,j,k) - fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
-
- enddo
- enddo
- enddo
-
- endif ! if (ispec_is_inner(ispec) .eqv. phase_is_inner)
-
- enddo ! spectral element loop
-
-! adding source
- do isource = 1,NSOURCES
-
- if (ispec_is_inner(ispec_selected_source(isource)) .eqv. phase_is_inner) then
-
- if(USE_FORCE_POINT_SOURCE) then
-
-! add the source (only if this proc carries the source)
- if(myrank == islice_selected_source(isource)) then
-
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
- t0 = 1.2d0/f0
-
- if (it == 1 .and. myrank == 0) then
- print *,'using a source of dominant frequency ',f0
- print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- endif
-
- ! we use nu_source(:,3) here because we want a source normal to the surface.
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- !accel(:,iglob) = accel(:,iglob) + &
- ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
- ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
- accel(:,iglob) = accel(:,iglob) + &
- sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
- exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
-
- endif
- endif
-
- endif
-
- enddo
-
-end subroutine compute_forces
-
-
-!! DK DK subroutines adapted from Deville, Fischer and Mund, High-order methods
-!! DK DK for incompressible fluid flow, Cambridge University Press (2002),
-!! DK DK pages 386 and 389 and Figure 8.3.1
-
- subroutine mxm_m1_m2_5points(A,B1,B2,B3,C1,C2,C3)
-
- implicit none
-
- include "constants.h"
-
- real(kind=CUSTOM_REAL), dimension(m1,NGLLX) :: A
- real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1,B2,B3
- real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1,C2,C3
-
- integer :: i,j
-
- do j=1,m2
- do i=1,m1
-
- C1(i,j) = A(i,1)*B1(1,j) + &
- A(i,2)*B1(2,j) + &
- A(i,3)*B1(3,j) + &
- A(i,4)*B1(4,j) + &
- A(i,5)*B1(5,j)
-
- C2(i,j) = A(i,1)*B2(1,j) + &
- A(i,2)*B2(2,j) + &
- A(i,3)*B2(3,j) + &
- A(i,4)*B2(4,j) + &
- A(i,5)*B2(5,j)
-
- C3(i,j) = A(i,1)*B3(1,j) + &
- A(i,2)*B3(2,j) + &
- A(i,3)*B3(3,j) + &
- A(i,4)*B3(4,j) + &
- A(i,5)*B3(5,j)
-
- enddo
- enddo
-
- end subroutine mxm_m1_m2_5points
-
-!---------
-
- subroutine mxm_m1_m1_5points(A1,A2,A3,B,C1,C2,C3)
-
- implicit none
-
- include "constants.h"
-
- real(kind=CUSTOM_REAL), dimension(m1,NGLLX) :: A1,A2,A3
- real(kind=CUSTOM_REAL), dimension(NGLLX,m1) :: B
- real(kind=CUSTOM_REAL), dimension(m1,m1) :: C1,C2,C3
-
- integer :: i,j
-
- do j=1,m1
- do i=1,m1
-
- C1(i,j) = A1(i,1)*B(1,j) + &
- A1(i,2)*B(2,j) + &
- A1(i,3)*B(3,j) + &
- A1(i,4)*B(4,j) + &
- A1(i,5)*B(5,j)
-
- C2(i,j) = A2(i,1)*B(1,j) + &
- A2(i,2)*B(2,j) + &
- A2(i,3)*B(3,j) + &
- A2(i,4)*B(4,j) + &
- A2(i,5)*B(5,j)
-
- C3(i,j) = A3(i,1)*B(1,j) + &
- A3(i,2)*B(2,j) + &
- A3(i,3)*B(3,j) + &
- A3(i,4)*B(4,j) + &
- A3(i,5)*B(5,j)
-
- enddo
- enddo
-
- end subroutine mxm_m1_m1_5points
-
-!---------
-
- subroutine mxm_m2_m1_5points(A1,A2,A3,B,C1,C2,C3)
-
- implicit none
-
- include "constants.h"
-
- real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1,A2,A3
- real(kind=CUSTOM_REAL), dimension(NGLLX,m1) :: B
- real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1,C2,C3
-
- integer :: i,j
-
- do j=1,m1
- do i=1,m2
-
- C1(i,j) = A1(i,1)*B(1,j) + &
- A1(i,2)*B(2,j) + &
- A1(i,3)*B(3,j) + &
- A1(i,4)*B(4,j) + &
- A1(i,5)*B(5,j)
-
- C2(i,j) = A2(i,1)*B(1,j) + &
- A2(i,2)*B(2,j) + &
- A2(i,3)*B(3,j) + &
- A2(i,4)*B(4,j) + &
- A2(i,5)*B(5,j)
-
- C3(i,j) = A3(i,1)*B(1,j) + &
- A3(i,2)*B(2,j) + &
- A3(i,3)*B(3,j) + &
- A3(i,4)*B(4,j) + &
- A3(i,5)*B(5,j)
-
- enddo
- enddo
-
- end subroutine mxm_m2_m1_5points
-
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -23,7 +23,7 @@
!
!=====================================================================
-subroutine compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+subroutine compute_forces_no_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
@@ -296,5 +296,5 @@
enddo
-end subroutine compute_forces
+end subroutine compute_forces_no_Deville
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 (from rev 15041, seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces.f90)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -0,0 +1,382 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+subroutine compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ kappastore,mustore,jacobian
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT
+ 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), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
+ newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! source
+ integer :: NSOURCES,myrank,it
+ integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
+ double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
+ double precision, dimension(3,3,NSOURCES) :: nu_source
+ double precision, dimension(NSOURCES) :: hdur
+ double precision :: dt
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+
+ integer ispec,iglob
+ integer i,j,k
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal
+
+ integer :: isource
+ double precision :: t0,f0
+
+ do ispec = 1,NSPEC_AB
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob)
+ enddo
+ enddo
+ enddo
+
+!! DK DK subroutines adapted from Deville, Fischer and Mund, High-order methods
+!! DK DK for incompressible fluid flow, Cambridge University Press (2002),
+!! DK DK pages 386 and 389 and Figure 8.3.1
+ call mxm_m1_m2_5points(hprime_xx,dummyx_loc,dummyy_loc,dummyz_loc,tempx1,tempy1,tempz1)
+
+ do k = 1,NGLLX
+ call mxm_m1_m1_5points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
+ hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
+ enddo
+
+ call mxm_m2_m1_5points(dummyx_loc,dummyy_loc,dummyz_loc,hprime_xxT,tempx3,tempy3,tempz3)
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ kappal = kappastore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
+
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+! form dot product with test vector, symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ enddo
+ enddo
+ enddo
+
+!! DK DK subroutines adapted from Deville, Fischer and Mund, High-order methods
+!! DK DK for incompressible fluid flow, Cambridge University Press (2002),
+!! DK DK pages 386 and 389 and Figure 8.3.1
+ call mxm_m1_m2_5points(hprimewgll_xxT,tempx1,tempy1,tempz1,newtempx1,newtempy1,newtempz1)
+
+ do k = 1,NGLLX
+ call mxm_m1_m1_5points(tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k), &
+ hprimewgll_xx,newtempx2(1,1,k),newtempy2(1,1,k),newtempz2(1,1,k))
+ enddo
+
+ call mxm_m2_m1_5points(tempx3,tempy3,tempz3,hprimewgll_xx,newtempx3,newtempy3,newtempz3)
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+! sum contributions from each element to the global mesh using indirect addressing
+ iglob = ibool(i,j,k,ispec)
+ accel(1,iglob) = accel(1,iglob) - fac1*newtempx1(i,j,k) - fac2*newtempx2(i,j,k) - fac3*newtempx3(i,j,k)
+ accel(2,iglob) = accel(2,iglob) - fac1*newtempy1(i,j,k) - fac2*newtempy2(i,j,k) - fac3*newtempy3(i,j,k)
+ accel(3,iglob) = accel(3,iglob) - fac1*newtempz1(i,j,k) - fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
+
+ enddo
+ enddo
+ enddo
+
+ endif ! if (ispec_is_inner(ispec) .eqv. phase_is_inner)
+
+ enddo ! spectral element loop
+
+! adding source
+ do isource = 1,NSOURCES
+
+ if (ispec_is_inner(ispec_selected_source(isource)) .eqv. phase_is_inner) then
+
+ if(USE_FORCE_POINT_SOURCE) then
+
+! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ t0 = 1.2d0/f0
+
+ if (it == 1 .and. myrank == 0) then
+ print *,'using a source of dominant frequency ',f0
+ print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ endif
+
+ ! we use nu_source(:,3) here because we want a source normal to the surface.
+ ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
+ !accel(:,iglob) = accel(:,iglob) + &
+ ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+ accel(:,iglob) = accel(:,iglob) + &
+ sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+
+ endif
+ endif
+
+ endif
+
+ enddo
+
+end subroutine compute_forces_with_Deville
+
+
+!! DK DK subroutines adapted from Deville, Fischer and Mund, High-order methods
+!! DK DK for incompressible fluid flow, Cambridge University Press (2002),
+!! DK DK pages 386 and 389 and Figure 8.3.1
+
+ subroutine mxm_m1_m2_5points(A,B1,B2,B3,C1,C2,C3)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL), dimension(m1,NGLLX) :: A
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1,B2,B3
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1,C2,C3
+
+ integer :: i,j
+
+ do j=1,m2
+ do i=1,m1
+
+ C1(i,j) = A(i,1)*B1(1,j) + &
+ A(i,2)*B1(2,j) + &
+ A(i,3)*B1(3,j) + &
+ A(i,4)*B1(4,j) + &
+ A(i,5)*B1(5,j)
+
+ C2(i,j) = A(i,1)*B2(1,j) + &
+ A(i,2)*B2(2,j) + &
+ A(i,3)*B2(3,j) + &
+ A(i,4)*B2(4,j) + &
+ A(i,5)*B2(5,j)
+
+ C3(i,j) = A(i,1)*B3(1,j) + &
+ A(i,2)*B3(2,j) + &
+ A(i,3)*B3(3,j) + &
+ A(i,4)*B3(4,j) + &
+ A(i,5)*B3(5,j)
+
+ enddo
+ enddo
+
+ end subroutine mxm_m1_m2_5points
+
+!---------
+
+ subroutine mxm_m1_m1_5points(A1,A2,A3,B,C1,C2,C3)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL), dimension(m1,NGLLX) :: A1,A2,A3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m1) :: B
+ real(kind=CUSTOM_REAL), dimension(m1,m1) :: C1,C2,C3
+
+ integer :: i,j
+
+ do j=1,m1
+ do i=1,m1
+
+ C1(i,j) = A1(i,1)*B(1,j) + &
+ A1(i,2)*B(2,j) + &
+ A1(i,3)*B(3,j) + &
+ A1(i,4)*B(4,j) + &
+ A1(i,5)*B(5,j)
+
+ C2(i,j) = A2(i,1)*B(1,j) + &
+ A2(i,2)*B(2,j) + &
+ A2(i,3)*B(3,j) + &
+ A2(i,4)*B(4,j) + &
+ A2(i,5)*B(5,j)
+
+ C3(i,j) = A3(i,1)*B(1,j) + &
+ A3(i,2)*B(2,j) + &
+ A3(i,3)*B(3,j) + &
+ A3(i,4)*B(4,j) + &
+ A3(i,5)*B(5,j)
+
+ enddo
+ enddo
+
+ end subroutine mxm_m1_m1_5points
+
+!---------
+
+ subroutine mxm_m2_m1_5points(A1,A2,A3,B,C1,C2,C3)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1,A2,A3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m1) :: B
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1,C2,C3
+
+ integer :: i,j
+
+ do j=1,m1
+ do i=1,m2
+
+ C1(i,j) = A1(i,1)*B(1,j) + &
+ A1(i,2)*B(2,j) + &
+ A1(i,3)*B(3,j) + &
+ A1(i,4)*B(4,j) + &
+ A1(i,5)*B(5,j)
+
+ C2(i,j) = A2(i,1)*B(1,j) + &
+ A2(i,2)*B(2,j) + &
+ A2(i,3)*B(3,j) + &
+ A2(i,4)*B(4,j) + &
+ A2(i,5)*B(5,j)
+
+ C3(i,j) = A3(i,1)*B(1,j) + &
+ A3(i,2)*B(2,j) + &
+ A3(i,3)*B(3,j) + &
+ A3(i,4)*B(4,j) + &
+ A3(i,5)*B(5,j)
+
+ enddo
+ enddo
+
+ end subroutine mxm_m2_m1_5points
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/configure
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/configure 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/configure 2009-05-25 16:12:48 UTC (rev 15043)
@@ -4988,7 +4988,7 @@
if test x"$LOCAL_PATH_IS_ALSO_GLOBAL" = x; then
- LOCAL_PATH_IS_ALSO_GLOBAL=false
+ LOCAL_PATH_IS_ALSO_GLOBAL=true
fi
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac 2009-05-25 16:12:48 UTC (rev 15043)
@@ -100,9 +100,9 @@
MPICC=mpicc
fi
-AC_ARG_VAR(LOCAL_PATH_IS_ALSO_GLOBAL, [files on a local path on each node are also seen as global with same path @<:@default=false@:>@])
+AC_ARG_VAR(LOCAL_PATH_IS_ALSO_GLOBAL, [files on a local path on each node are also seen as global with same path @<:@default=true@:>@])
if test x"$LOCAL_PATH_IS_ALSO_GLOBAL" = x; then
- LOCAL_PATH_IS_ALSO_GLOBAL=false
+ LOCAL_PATH_IS_ALSO_GLOBAL=true
fi
AC_ARG_VAR(PYTHON, [Python interpreter])
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2009-05-25 16:12:48 UTC (rev 15043)
@@ -44,6 +44,9 @@
! if running on a Beowulf, also modify name of nodes in filter_machine_file.f90
logical, parameter :: LOCAL_PATH_IS_ALSO_GLOBAL = . at LOCAL_PATH_IS_ALSO_GLOBAL@.
+! use inlined products of Deville et al. (2002) to speedup the calculations to compute internal forces
+ logical, parameter :: USE_DEVILLE_PRODUCTS = .true.
+
! apply heuristic rule to modify doubling regions to balance angles
logical, parameter :: APPLY_HEURISTIC_RULE = .true.
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_header_file.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_header_file.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -96,8 +96,11 @@
print *
print *,'edit file OUTPUT_FILES/values_from_mesher.h to see some statistics about the mesh'
print *
- print *,'on NEC SX, make sure "loopcnt=" parameter'
- print *,'in Makefile is greater than max vector length = ',NGLOB_AB
+!! DK DK May 2009: removed this because now each slice of a CUBIT + SCOTCH mesh
+!! DK DK May 2009: has a different number of spectral elements and therefore the
+!! DK DK May 2009: value below should be the max() for all the slices
+! print *,'on NEC SX, make sure "loopcnt=" parameter'
+! print *,'in Makefile is greater than max vector length = ',NGLOB_AB
print *
print *,'done'
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -1255,6 +1255,7 @@
integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
integer :: inumber
+
! **************
! create the name for the database of the current slide and region
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -59,15 +59,29 @@
write(IOUT,*) '! mesh statistics:'
write(IOUT,*) '! ---------------'
write(IOUT,*) '!'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK these statistics are now INCORRECT'
+ write(IOUT,*) '! DK DK because the CUBIT + SCOTCH mesh has'
+ write(IOUT,*) '! DK DK a different number of mesh elements and points in each slice'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '! DK DK'
+ write(IOUT,*) '!'
write(IOUT,*) '! number of processors = ',NPROC
write(IOUT,*) '!'
write(IOUT,*) '! number of ES nodes = ',real(NPROC)/8.
write(IOUT,*) '! percentage of total 640 ES nodes = ',100.*(real(NPROC)/8.)/640.,' %'
write(IOUT,*) '! total memory available on these ES nodes (Gb) = ',16.*real(NPROC)/8.
- write(IOUT,*) 'integer, parameter :: NPROC_VAL = ',NPROC
- write(IOUT,*) 'integer, parameter :: NPROC_XI_VAL = ', NPROC_XI
- write(IOUT,*) 'integer, parameter :: NPROC_ETA_VAL = ', NPROC_ETA
+! write(IOUT,*) 'integer, parameter :: NPROC_VAL = ',NPROC
+! write(IOUT,*) 'integer, parameter :: NPROC_XI_VAL = ', NPROC_XI
+! write(IOUT,*) 'integer, parameter :: NPROC_ETA_VAL = ', NPROC_ETA
write(IOUT,*) '!'
write(IOUT,*) '! max points per processor = max vector length = ',NGLOB_AB
@@ -104,16 +118,12 @@
write(IOUT,*) '! average distance between points along X in m = ',sngl(UTM_X_MAX-UTM_X_MIN)/real(NEX_XI*(NGLLX-1))
write(IOUT,*) '! average distance between points along Y in m = ',sngl(UTM_Y_MAX-UTM_Y_MIN)/real(NEX_ETA*(NGLLY-1))
write(IOUT,*) '!'
- write(IOUT,*)
- write(IOUT,*) 'integer, parameter :: NSPEC_AB_VAL = ',NSPEC_AB
- write(IOUT,*) 'integer, parameter :: NGLOB_AB_VAL = ',NGLOB_AB
- write(IOUT,*)
- write(IOUT,*) '!'
write(IOUT,*) '! number of time steps = ',NSTEP
write(IOUT,*) '!'
! if attenuation is off, set dummy size of arrays to one
if(ATTENUATION) then
+ stop 'ATTENUATION not supported yet in the CUBIT + SCOTCH version because of arrays of constant size defined'
write(IOUT,*) 'integer, parameter :: NSPEC_ATTENUATION = ', NSPEC_AB
write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .true.'
else
@@ -124,6 +134,7 @@
! anisotropy
if(ANISOTROPY) then
+ stop 'ANISOTROPY not supported yet in the CUBIT + SCOTCH version because of arrays of constant size defined'
write(IOUT,*) 'integer, parameter :: NSPEC_ANISO = ',NSPEC_AB
write(IOUT,*) 'logical, parameter :: ANISOTROPY_VAL = .true.'
else
@@ -133,41 +144,45 @@
write(IOUT,*)
+!! DK DK May 2009: removed all the things that are not supported in the CUBIT + SCOTCH version yet
+!! DK DK May 2009: removed all the things that are not supported in the CUBIT + SCOTCH version yet
+!! DK DK May 2009: removed all the things that are not supported in the CUBIT + SCOTCH version yet
+
! strain/attenuation
if (ATTENUATION .and. SIMULATION_TYPE == 3) then
- write(IOUT,*) 'integer, parameter :: NSPEC_ATT_AND_KERNEL = ', NSPEC_AB
+! write(IOUT,*) 'integer, parameter :: NSPEC_ATT_AND_KERNEL = ', NSPEC_AB
else
- write(IOUT,*) 'integer, parameter :: NSPEC_ATT_AND_KERNEL = ', 1
+! write(IOUT,*) 'integer, parameter :: NSPEC_ATT_AND_KERNEL = ', 1
endif
! adjoint
if (SIMULATION_TYPE == 3) then
- write(IOUT,*) 'integer, parameter :: NSPEC_ADJOINT = ', NSPEC_AB
- write(IOUT,*) 'integer, parameter :: NGLOB_ADJOINT = ', NGLOB_AB
+! write(IOUT,*) 'integer, parameter :: NSPEC_ADJOINT = ', NSPEC_AB
+! write(IOUT,*) 'integer, parameter :: NGLOB_ADJOINT = ', NGLOB_AB
else
- write(IOUT,*) 'integer, parameter :: NSPEC_ADJOINT = ', 1
- write(IOUT,*) 'integer, parameter :: NGLOB_ADJOINT = ', 1
+! write(IOUT,*) 'integer, parameter :: NSPEC_ADJOINT = ', 1
+! write(IOUT,*) 'integer, parameter :: NGLOB_ADJOINT = ', 1
endif
write(IOUT,*)
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_VAL = ', NSPEC2DMAX_XMIN_XMAX
- write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_VAL = ', NSPEC2DMAX_YMIN_YMAX
- write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_VAL = ', NSPEC2D_BOTTOM
- write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_VAL = ', NSPEC2D_TOP
- write(IOUT,*) 'integer, parameter :: NPOIN2DMAX_XMIN_XMAX_VAL = ', NPOIN2DMAX_XMIN_XMAX
- write(IOUT,*) 'integer, parameter :: NPOIN2DMAX_YMIN_YMAX_VAL = ', NPOIN2DMAX_YMIN_YMAX
- write(IOUT,*) 'integer, parameter :: NPOIN2DMAX_XY_VAL = ', NPOIN2DMAX_XY
+! write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_XMIN_XMAX_VAL = ', NSPEC2DMAX_XMIN_XMAX
+! write(IOUT,*) 'integer, parameter :: NSPEC2DMAX_YMIN_YMAX_VAL = ', NSPEC2DMAX_YMIN_YMAX
+! write(IOUT,*) 'integer, parameter :: NSPEC2D_BOTTOM_VAL = ', NSPEC2D_BOTTOM
+! write(IOUT,*) 'integer, parameter :: NSPEC2D_TOP_VAL = ', NSPEC2D_TOP
+! write(IOUT,*) 'integer, parameter :: NPOIN2DMAX_XMIN_XMAX_VAL = ', NPOIN2DMAX_XMIN_XMAX
+! write(IOUT,*) 'integer, parameter :: NPOIN2DMAX_YMIN_YMAX_VAL = ', NPOIN2DMAX_YMIN_YMAX
+! write(IOUT,*) 'integer, parameter :: NPOIN2DMAX_XY_VAL = ', NPOIN2DMAX_XY
write(IOUT,*)
! Moho boundary
if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
- write(IOUT,*) 'integer, parameter :: NSPEC2D_MOHO_BOUN = ', NSPEC2D_BOTTOM
- write(IOUT,*) 'integer, parameter :: NSPEC_BOUN = ', NSPEC_AB
+! write(IOUT,*) 'integer, parameter :: NSPEC2D_MOHO_BOUN = ', NSPEC2D_BOTTOM
+! write(IOUT,*) 'integer, parameter :: NSPEC_BOUN = ', NSPEC_AB
else
- write(IOUT,*) 'integer, parameter :: NSPEC2D_MOHO_BOUN = ', 1
- write(IOUT,*) 'integer, parameter :: NSPEC_BOUN = ', 1
+! write(IOUT,*) 'integer, parameter :: NSPEC2D_MOHO_BOUN = ', 1
+! write(IOUT,*) 'integer, parameter :: NSPEC_BOUN = ', 1
endif
close(IOUT)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90 2009-05-25 15:55:42 UTC (rev 15042)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90 2009-05-25 16:12:48 UTC (rev 15043)
@@ -34,6 +34,9 @@
! include values created by the mesher
include "OUTPUT_FILES/values_from_mesher.h"
+! standard include of the MPI library
+ include 'mpif.h'
+
!=============================================================================!
! !
! specfem3D is a 3-D spectral-element solver for a local or regional model. !
@@ -145,10 +148,11 @@
! ADJOINT
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: b_alphaval, b_betaval, b_gammaval
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) :: &
- b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL) :: b_epsilondev_xx, &
- b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
+!! DK DK array not created yet for CUBIT
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) :: &
+! b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL) :: b_epsilondev_xx, &
+! b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
! ADJOINT
integer NPOIN2DMAX_XY
@@ -159,16 +163,19 @@
character(len=100) topo_file
integer, dimension(:,:), allocatable :: itopo_bathy
- integer, dimension(NSPEC2D_TOP_VAL) :: ibelm_top
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_VAL) :: normal_top
+!! DK DK array not created yet for CUBIT
+! integer, dimension(NSPEC2D_TOP_VAL) :: ibelm_top
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_VAL) :: normal_top
+!! DK DK array not created yet for CUBIT
! Moho mesh
- integer,dimension(NSPEC2D_MOHO_BOUN) :: ibelm_moho_top, ibelm_moho_bot
- real(CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_MOHO_BOUN) :: normal_moho
+! integer,dimension(NSPEC2D_MOHO_BOUN) :: ibelm_moho_top, ibelm_moho_bot
+! real(CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_MOHO_BOUN) :: normal_moho
integer :: nspec2D_moho
+!! DK DK array not created yet for CUBIT
! buffers for send and receive between faces of the slices and the chunks
- real(kind=CUSTOM_REAL), dimension(NDIM,NPOIN2DMAX_XY_VAL) :: buffer_send_faces_vector,buffer_received_faces_vector
+! real(kind=CUSTOM_REAL), dimension(NDIM,NPOIN2DMAX_XY_VAL) :: buffer_send_faces_vector,buffer_received_faces_vector
! -----------------
@@ -217,9 +224,10 @@
! ADJOINT
real(kind=CUSTOM_REAL) b_additional_term,b_force_normal_comp
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT) :: b_displ, b_veloc, b_accel
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: rho_kl, mu_kl, kappa_kl, &
- rhop_kl, beta_kl, alpha_kl
+!! DK DK array not created yet for CUBIT
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT) :: b_displ, b_veloc, b_accel
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: rho_kl, mu_kl, kappa_kl, &
+! rhop_kl, beta_kl, alpha_kl
real(kind=CUSTOM_REAL) b_deltat, b_deltatover2, b_deltatsqover2
! ADJOINT
@@ -229,8 +237,9 @@
! Moho kernel
integer ispec2D_moho_top, ispec2D_moho_bot, k_top, k_bot, ispec_top, ispec_bot, iglob_top, iglob_bot
- real(kind=CUSTOM_REAL), dimension(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO_BOUN) :: dsdx_top, dsdx_bot, b_dsdx_top, b_dsdx_bot
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_MOHO_BOUN) :: moho_kl
+!! DK DK array not created yet for CUBIT
+! real(kind=CUSTOM_REAL), dimension(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO_BOUN) :: dsdx_top, dsdx_bot, b_dsdx_top, b_dsdx_bot
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_MOHO_BOUN) :: moho_kl
real(kind=CUSTOM_REAL) :: kernel_moho_top, kernel_moho_bot
! --------
@@ -300,11 +309,11 @@
double precision, dimension(:,:), allocatable :: hxir_store,hetar_store,hgammar_store
! 2-D addressing and buffers for summation between slices
- integer, dimension(NPOIN2DMAX_XMIN_XMAX_VAL) :: iboolleft_xi,iboolright_xi
- integer, dimension(NPOIN2DMAX_YMIN_YMAX_VAL) :: iboolleft_eta,iboolright_eta
+! integer, dimension(NPOIN2DMAX_XMIN_XMAX_VAL) :: iboolleft_xi,iboolright_xi
+! integer, dimension(NPOIN2DMAX_YMIN_YMAX_VAL) :: iboolleft_eta,iboolright_eta
! for addressing of the slices
- integer, dimension(0:NPROC_XI_VAL-1,0:NPROC_ETA_VAL) :: addressing
+! integer, dimension(0:NPROC_XI_VAL-1,0:NPROC_ETA_VAL) :: addressing
! proc numbers for MPI
integer myrank,sizeprocs
@@ -428,6 +437,12 @@
!!$ real(kind=CUSTOM_REAL) :: weight, jacobianl
!!!! NL NL REGOLITH
+!! DK DK May 2009: added this to print the minimum and maximum number of elements
+!! DK DK May 2009: and points in the CUBIT + SCOTCH mesh
+ integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum
+ integer :: NGLOB_AB_global_min,NGLOB_AB_global_max
+ integer :: ier
+
! ************** PROGRAM STARTS HERE **************
! sizeprocs returns number of processes started
@@ -1543,18 +1558,19 @@
! put negligible initial value to avoid very slow underflow trapping
if(FIX_UNDERFLOW_PROBLEM) displ(:,:) = VERYSMALLVAL
- if (SIMULATION_TYPE == 3) then ! kernel calculation, read in last frame
+!! DK DK array not created yet for CUBIT
+! if (SIMULATION_TYPE == 3) then ! kernel calculation, read in last frame
- open(unit=27,file=trim(prname)//'save_forward_arrays.bin',status='old',action='read',form='unformatted')
- read(27) b_displ
- read(27) b_veloc
- read(27) b_accel
+! open(unit=27,file=trim(prname)//'save_forward_arrays.bin',status='old',action='read',form='unformatted')
+! read(27) b_displ
+! read(27) b_veloc
+! read(27) b_accel
- rho_kl(:,:,:,:) = 0._CUSTOM_REAL
- mu_kl(:,:,:,:) = 0._CUSTOM_REAL
- kappa_kl(:,:,:,:) = 0._CUSTOM_REAL
+! rho_kl(:,:,:,:) = 0._CUSTOM_REAL
+! mu_kl(:,:,:,:) = 0._CUSTOM_REAL
+! kappa_kl(:,:,:,:) = 0._CUSTOM_REAL
- endif
+! endif
! allocate files to save movies and shaking map
if(MOVIE_SURFACE .or. CREATE_SHAKEMAP) then
@@ -1667,18 +1683,19 @@
R_yz(:,:,:,:,:) = VERYSMALLVAL
endif
- if (SIMULATION_TYPE == 3) then
- read(27) b_R_xx
- read(27) b_R_yy
- read(27) b_R_xy
- read(27) b_R_xz
- read(27) b_R_yz
- read(27) b_epsilondev_xx
- read(27) b_epsilondev_yy
- read(27) b_epsilondev_xy
- read(27) b_epsilondev_xz
- read(27) b_epsilondev_yz
- endif
+!! DK DK array not created yet for CUBIT
+! if (SIMULATION_TYPE == 3) then
+! read(27) b_R_xx
+! read(27) b_R_yy
+! read(27) b_R_xy
+! read(27) b_R_xz
+! read(27) b_R_yz
+! read(27) b_epsilondev_xx
+! read(27) b_epsilondev_yy
+! read(27) b_epsilondev_xy
+! read(27) b_epsilondev_xz
+! read(27) b_epsilondev_yz
+! endif
endif
close(27)
@@ -1691,6 +1708,29 @@
k_bot = NGLLZ
endif
+!! DK DK May 2009: added this to print the minimum and maximum number of elements
+!! DK DK May 2009: and points in the CUBIT + SCOTCH mesh
+ call MPI_REDUCE(NSPEC_AB,NSPEC_AB_global_min,1,MPI_INTEGER,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call MPI_REDUCE(NSPEC_AB,NSPEC_AB_global_max,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call MPI_REDUCE(NSPEC_AB,NSPEC_AB_global_sum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ call MPI_REDUCE(NGLOB_AB,NGLOB_AB_global_min,1,MPI_INTEGER,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call MPI_REDUCE(NGLOB_AB,NGLOB_AB_global_max,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'minimum and maximum number of elements'
+ write(IMAIN,*) 'and points in the CUBIT + SCOTCH mesh:'
+ write(IMAIN,*)
+ write(IMAIN,*) 'NSPEC_AB_global_min = ',NSPEC_AB_global_min
+ write(IMAIN,*) 'NSPEC_AB_global_max = ',NSPEC_AB_global_max
+ write(IMAIN,*) 'NSPEC_AB_global_mean = ',NSPEC_AB_global_sum / float(sizeprocs)
+ write(IMAIN,*)
+ write(IMAIN,*) 'NGLOB_AB_global_min = ',NGLOB_AB_global_min
+ write(IMAIN,*) 'NGLOB_AB_global_max = ',NGLOB_AB_global_max
+ write(IMAIN,*)
+ endif
+
!
! s t a r t t i m e i t e r a t i o n s
!
@@ -1732,10 +1772,12 @@
! compute the maximum of the maxima for all the slices using an MPI reduction
call max_all_cr(Usolidnorm,Usolidnorm_all)
- if (SIMULATION_TYPE == 3) then
- b_Usolidnorm = maxval(sqrt(b_displ(1,:)**2 + b_displ(2,:)**2 + b_displ(3,:)**2))
- call max_all_cr(b_Usolidnorm,b_Usolidnorm_all)
- endif
+!! DK DK array not created yet for CUBIT
+! if (SIMULATION_TYPE == 3) then
+! b_Usolidnorm = maxval(sqrt(b_displ(1,:)**2 + b_displ(2,:)**2 + b_displ(3,:)**2))
+! call max_all_cr(b_Usolidnorm,b_Usolidnorm_all)
+! endif
+
if(myrank == 0) then
write(IMAIN,*) 'Time step # ',it
@@ -1784,11 +1826,12 @@
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
accel(:,:) = 0._CUSTOM_REAL
- if (SIMULATION_TYPE == 3) then
- b_displ(:,:) = b_displ(:,:) + b_deltat*b_veloc(:,:) + b_deltatsqover2*b_accel(:,:)
- b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
- b_accel(:,:) = 0._CUSTOM_REAL
- endif
+!! DK DK array not created yet for CUBIT
+! if (SIMULATION_TYPE == 3) then
+! b_displ(:,:) = b_displ(:,:) + b_deltat*b_veloc(:,:) + b_deltatsqover2*b_accel(:,:)
+! b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
+! b_accel(:,:) = 0._CUSTOM_REAL
+! endif
if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
ispec2D_moho_top = 0
@@ -1796,10 +1839,17 @@
endif
! assemble all the contributions between slices using MPI
- call compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ if(USE_DEVILLE_PRODUCTS) then
+ call compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.false., &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ else
+ call compute_forces_no_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.false., &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ endif
call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_AB,accel, &
buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
@@ -1807,41 +1857,56 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
- call compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ if(USE_DEVILLE_PRODUCTS) then
+ call compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.true., &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ else
+ call compute_forces_no_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.true., &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ endif
call assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_AB,accel, &
buffer_recv_vector_ext_mesh,ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
- if (SIMULATION_TYPE == 3) call assemble_MPI_vector(b_accel,iproc_xi,iproc_eta,addressing, &
- iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
- buffer_send_faces_vector,buffer_received_faces_vector,npoin2D_xi,npoin2D_eta, &
- NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
+!! DK DK May 2009: removed this because now each slice of a CUBIT + SCOTCH mesh
+!! DK DK May 2009: has a different number of spectral elements and therefore
+!! DK DK May 2009: only the general non-blocking MPI routines assemble_MPI_vector_ext_mesh_s
+!! DK DK May 2009: and assemble_MPI_vector_ext_mesh_w above can be used.
+!! DK DK May 2009: For adjoint runs below (SIMULATION_TYPE == 3) they should be used as well.
+! if (SIMULATION_TYPE == 3) call assemble_MPI_vector(b_accel,iproc_xi,iproc_eta,addressing, &
+! iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
+! buffer_send_faces_vector,buffer_received_faces_vector,npoin2D_xi,npoin2D_eta, &
+! NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
! multiply by the inverse of the mass matrix
accel(1,:) = accel(1,:)*rmass(:)
accel(2,:) = accel(2,:)*rmass(:)
accel(3,:) = accel(3,:)*rmass(:)
- if (SIMULATION_TYPE == 3) then
- b_accel(1,:) = b_accel(1,:)*rmass(:)
- b_accel(2,:) = b_accel(2,:)*rmass(:)
- b_accel(3,:) = b_accel(3,:)*rmass(:)
- endif
+!! DK DK array not created yet for CUBIT
+! if (SIMULATION_TYPE == 3) then
+! b_accel(1,:) = b_accel(1,:)*rmass(:)
+! b_accel(2,:) = b_accel(2,:)*rmass(:)
+! b_accel(3,:) = b_accel(3,:)*rmass(:)
+! endif
if(OCEANS) then
+ stop 'DK DK oceans have been removed for now because we need a flag to detect the surface elements'
+
! initialize the updates
updated_dof_ocean_load(:) = .false.
! for surface elements exactly at the top of the model (ocean bottom)
do ispec2D = 1,NSPEC2D_TOP
- ispec = ibelm_top(ispec2D)
+!! DK DK array not created yet for CUBIT ispec = ibelm_top(ispec2D)
! only for DOFs exactly at the top of the model (ocean bottom)
k = NGLLZ
@@ -1856,9 +1921,9 @@
if(.not. updated_dof_ocean_load(iglob)) then
! get normal
- nx = normal_top(1,i,j,ispec2D)
- ny = normal_top(2,i,j,ispec2D)
- nz = normal_top(3,i,j,ispec2D)
+!! DK DK array not created yet for CUBIT nx = normal_top(1,i,j,ispec2D)
+!! DK DK array not created yet for CUBIT ny = normal_top(2,i,j,ispec2D)
+!! DK DK array not created yet for CUBIT nz = normal_top(3,i,j,ispec2D)
! make updated component of right-hand side
! we divide by rmass() which is 1 / M
@@ -1873,14 +1938,16 @@
accel(3,iglob) = accel(3,iglob) + additional_term * nz
if (SIMULATION_TYPE == 3) then
- b_force_normal_comp = (b_accel(1,iglob)*nx + &
- b_accel(2,iglob)*ny + b_accel(3,iglob)*nz) / rmass(iglob)
+!! DK DK array not created yet for CUBIT
+! b_force_normal_comp = (b_accel(1,iglob)*nx + &
+! b_accel(2,iglob)*ny + b_accel(3,iglob)*nz) / rmass(iglob)
b_additional_term = (rmass_ocean_load(iglob) - rmass(iglob)) * b_force_normal_comp
- b_accel(1,iglob) = b_accel(1,iglob) + b_additional_term * nx
- b_accel(2,iglob) = b_accel(2,iglob) + b_additional_term * ny
- b_accel(3,iglob) = b_accel(3,iglob) + b_additional_term * nz
+!! DK DK array not created yet for CUBIT
+! b_accel(1,iglob) = b_accel(1,iglob) + b_additional_term * nx
+! b_accel(2,iglob) = b_accel(2,iglob) + b_additional_term * ny
+! b_accel(3,iglob) = b_accel(3,iglob) + b_additional_term * nz
endif
! done with this point
@@ -1895,7 +1962,8 @@
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
- if (SIMULATION_TYPE == 3) b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
+!! DK DK array not created yet for CUBIT
+! if (SIMULATION_TYPE == 3) b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
! write the seismograms with time shift
if (nrec_local > 0) then
@@ -2019,15 +2087,16 @@
hlagrange = hxir_store(irec_local,i)*hetar_store(irec_local,j)*hgammar_store(irec_local,k)
- dxd = dxd + dble(b_displ(1,iglob))*hlagrange
- dyd = dyd + dble(b_displ(2,iglob))*hlagrange
- dzd = dzd + dble(b_displ(3,iglob))*hlagrange
- vxd = vxd + dble(b_veloc(1,iglob))*hlagrange
- vyd = vyd + dble(b_veloc(2,iglob))*hlagrange
- vzd = vzd + dble(b_veloc(3,iglob))*hlagrange
- axd = axd + dble(b_accel(1,iglob))*hlagrange
- ayd = ayd + dble(b_accel(2,iglob))*hlagrange
- azd = azd + dble(b_accel(3,iglob))*hlagrange
+!! DK DK array not created yet for CUBIT
+! dxd = dxd + dble(b_displ(1,iglob))*hlagrange
+! dyd = dyd + dble(b_displ(2,iglob))*hlagrange
+! dzd = dzd + dble(b_displ(3,iglob))*hlagrange
+! vxd = vxd + dble(b_veloc(1,iglob))*hlagrange
+! vyd = vyd + dble(b_veloc(2,iglob))*hlagrange
+! vzd = vzd + dble(b_veloc(3,iglob))*hlagrange
+! axd = axd + dble(b_accel(1,iglob))*hlagrange
+! ayd = ayd + dble(b_accel(2,iglob))*hlagrange
+! azd = azd + dble(b_accel(3,iglob))*hlagrange
enddo
enddo
enddo
@@ -2074,19 +2143,20 @@
if (SIMULATION_TYPE == 3 .and. mod(NSTEP-it,NSTEP_Q_SAVE) == 0) then
write(outputname,"('save_Q_arrays_',i6.6,'.bin')") NSTEP-it
open(unit=27,file=trim(prname_Q)//trim(outputname),status='old',action='read',form='unformatted')
- read(27) b_displ
- read(27) b_veloc
- read(27) b_accel
- read(27) b_R_xx
- read(27) b_R_yy
- read(27) b_R_xy
- read(27) b_R_xz
- read(27) b_R_yz
- read(27) b_epsilondev_xx
- read(27) b_epsilondev_yy
- read(27) b_epsilondev_xy
- read(27) b_epsilondev_xz
- read(27) b_epsilondev_yz
+!! DK DK array not created yet for CUBIT
+! read(27) b_displ
+! read(27) b_veloc
+! read(27) b_accel
+! read(27) b_R_xx
+! read(27) b_R_yy
+! read(27) b_R_xy
+! read(27) b_R_xz
+! read(27) b_R_yz
+! read(27) b_epsilondev_xx
+! read(27) b_epsilondev_yy
+! read(27) b_epsilondev_xy
+! read(27) b_epsilondev_xz
+! read(27) b_epsilondev_yz
close(27)
else if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. mod(it,NSTEP_Q_SAVE) == 0) then
write(outputname,"('save_Q_arrays_',i6.6,'.bin')") it
@@ -2108,49 +2178,6 @@
endif
endif
-! kernel calculations
- if (SIMULATION_TYPE == 3) then
- do ispec = 1, NSPEC_AB
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- rho_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) + deltat * dot_product(accel(:,iglob), b_displ(:,iglob))
- enddo
- enddo
- enddo
- enddo
-
- if (SAVE_MOHO_MESH) then
- do ispec2D = 1, nspec2D_moho
- ispec_top = ibelm_moho_top(ispec2D)
- ispec_bot = ibelm_moho_bot(ispec2D)
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob_top = ibool(i,j,k_top,ispec_top)
-
- call compute_boundary_kernel(kernel_moho_top, &
- mustore(i,j,k_top,ispec_top), kappastore(i,j,k_top,ispec_top), rho_vs(i,j,k_top,ispec_top), &
- accel(:,iglob_top),b_displ(:,iglob_top),dsdx_top(:,:,i,j,k_top,ispec2D), b_dsdx_top(:,:,i,j,k_top,ispec2D), &
- normal_moho(:,i,j,ispec2D))
-
- iglob_bot = ibool(i,j,k_bot,ispec_bot)
- ! iglob_top == iglob_bot!
-
- call compute_boundary_kernel(kernel_moho_bot, &
- mustore(i,j,k_bot,ispec_bot), kappastore(i,j,k_bot,ispec_bot), rho_vs(i,j,k_bot,ispec_bot), &
- accel(:,iglob_bot),b_displ(:,iglob_bot),dsdx_bot(:,:,i,j,k_bot,ispec2D), b_dsdx_bot(:,:,i,j,k_bot,ispec2D), &
- normal_moho(:,i,j,ispec2D))
-
- moho_kl(i,j,ispec2D) = moho_kl(i,j,ispec2D) + (kernel_moho_top - kernel_moho_bot) * deltat
-
- enddo
- enddo
- enddo
- endif
- endif
-
-
if (EXTERNAL_MESH_CREATE_SHAKEMAP) then
if (it == 1) then
@@ -2416,13 +2443,15 @@
! save MOVIE on the SURFACE
if(MOVIE_SURFACE .and. mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
+ stop 'DK DK MOVIE_SURFACE has been removed for now because we need a flag to detect the surface elements'
+
! get coordinates of surface mesh and surface displacement
ipoin = 0
k = NGLLZ
if (USE_HIGHRES_FOR_MOVIES) then
do ispec2D = 1,NSPEC2D_TOP
- ispec = ibelm_top(ispec2D)
+!! DK DK array not created yet for CUBIT ispec = ibelm_top(ispec2D)
do j = 1,NGLLY
do i = 1,NGLLX
ipoin = ipoin + 1
@@ -2444,7 +2473,7 @@
enddo ! ispec_top
else
do ispec2D = 1,NSPEC2D_TOP
- ispec = ibelm_top(ispec2D)
+!! DK DK array not created yet for CUBIT ispec = ibelm_top(ispec2D)
do iloc = 1, NGNOD2D
ipoin = ipoin + 1
iglob = ibool(iorderi(iloc),iorderj(iloc),k,ispec)
@@ -2491,6 +2520,8 @@
! compute SHAKING INTENSITY MAP
if(CREATE_SHAKEMAP) then
+ stop 'DK DK CREATE_SHAKEMAP has been removed for now because we need a flag to detect the surface elements'
+
ipoin = 0
k = NGLLZ
@@ -2498,7 +2529,7 @@
if(USE_HIGHRES_FOR_MOVIES) then
do ispec2D = 1,NSPEC2D_TOP
- ispec = ibelm_top(ispec2D)
+!! DK DK array not created yet for CUBIT ispec = ibelm_top(ispec2D)
! loop on all the points inside the element
do j = 1,NGLLY
@@ -2517,7 +2548,7 @@
else
do ispec2D = 1,NSPEC2D_TOP
- ispec = ibelm_top(ispec2D)
+!! DK DK array not created yet for CUBIT ispec = ibelm_top(ispec2D)
do iloc = 1, NGNOD2D
ipoin = ipoin + 1
iglob = ibool(iorderi(iloc),iorderj(iloc),k,ispec)
@@ -2700,49 +2731,8 @@
else if (SIMULATION_TYPE == 3) then
! rhop, beta, alpha kernels
- do ispec = 1, NSPEC_AB
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- rho_kl(i,j,k,ispec) = - rho_vs(i,j,k,ispec) **2 * rho_kl(i,j,k,ispec) / mustore(i,j,k,ispec)
- mu_kl(i,j,k,ispec) = - mustore(i,j,k,ispec) * mu_kl(i,j,k,ispec)
- kappa_kl(i,j,k,ispec) = - kappastore(i,j,k,ispec) * kappa_kl(i,j,k,ispec)
- rhop_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) + kappa_kl(i,j,k,ispec) + mu_kl(i,j,k,ispec)
- beta_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (mu_kl(i,j,k,ispec) - 4._CUSTOM_REAL * mustore(i,j,k,ispec) &
- / (3._CUSTOM_REAL * kappastore(i,j,k,ispec)) * kappa_kl(i,j,k,ispec))
- alpha_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (1._CUSTOM_REAL + &
- 4._CUSTOM_REAL * mustore(i,j,k,ispec)/ (3._CUSTOM_REAL * kappastore(i,j,k,ispec))) &
- * kappa_kl(i,j,k,ispec)
- enddo
- enddo
- enddo
- enddo
-
! save kernels to binary files
- open(unit=27,file=prname(1:len_trim(prname))//'rho_kernel.bin',status='unknown',form='unformatted')
- write(27) rho_kl
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'mu_kernel.bin',status='unknown',form='unformatted')
- write(27) mu_kl
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'kappa_kernel.bin',status='unknown',form='unformatted')
- write(27) kappa_kl
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'rhop_kernel.bin',status='unknown',form='unformatted')
- write(27) rhop_kl
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'beta_kernel.bin',status='unknown',form='unformatted')
- write(27) beta_kl
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'alpha_kernel.bin',status='unknown',form='unformatted')
- write(27) alpha_kl
- close(27)
- if (SAVE_MOHO_MESH) then
- open(unit=27,file=prname(1:len_trim(prname))//'moho_kernel.bin',status='unknown',form='unformatted')
- write(27) moho_kl
- close(27)
- endif
+!! DK DK removed kernels from here because not supported for CUBIT + SCOTCH yet
endif
More information about the CIG-COMMITS
mailing list