[cig-commits] r12671 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: . src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Aug 17 19:29:55 PDT 2008
Author: dkomati1
Date: 2008-08-17 19:29:55 -0700 (Sun, 17 Aug 2008)
New Revision: 12671
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90
Log:
fixed the "negative mass matrix term" bug, which was related to the fact that I had
removed the second pass in the mesher; in fact, my algorithm is designed for two
passes and cannot work with only one pass; I therefore put the second pass back
but I now assign model parameters (velocities and density) in the second pass only
to save CPU time in the case of 3D models that require a costly interpolation
process at grid points.
Also added -mcmodel=medium compiler option for cases that require more than 2 GB
of memory per core.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-08-18 00:39:47 UTC (rev 12670)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-08-18 02:29:55 UTC (rev 12671)
@@ -37,8 +37,8 @@
FC = ifort
MPIFC = mpif90
MPIFLAGS = -DUSE_MPI # -lmpi
-#FLAGS_NO_CHECK = -O1 -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe0 -ftz -traceback -ftrapuv
-FLAGS_NO_CHECK = -O3 -xP -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
+#FLAGS_NO_CHECK = -O1 -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe0 -ftz -traceback -ftrapuv # -mcmodel=medium -shared-intel
+FLAGS_NO_CHECK = -O3 -xP -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz # -mcmodel=medium -shared-intel
# we need the -no-heap-arrays flag to force the compiler to allocate memory on the stack
# instead of on the heap to minimize memory fragmentation
@@ -48,7 +48,7 @@
#FC = gfortran
#MPIFC = /opt/mpich2_gfortran/bin/mpif90
#MPIFLAGS = -DUSE_MPI
-#FLAGS_NO_CHECK = -std=f2003 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -fbounds-check
+#FLAGS_NO_CHECK = -std=f2003 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -fbounds-check # -mcmodel=medium
#
# Portland pgf90
@@ -56,7 +56,7 @@
#FC = pgf90
#MPIFC = mpif90
#MPIFLAGS = -DUSE_MPI
-#FLAGS_NO_CHECK = -fast -Mnobounds -Mrecursive -Minline -Mdclchk -Knoieee -fastsse -tp amd64e -Minform=warn -Ktrap=none
+#FLAGS_NO_CHECK = -fast -Mnobounds -Mrecursive -Minline -Mdclchk -Knoieee -fastsse -tp amd64e -Minform=warn -Ktrap=none # -mcmodel=medium
# we need the -Mrecursive flag to force the compiler to allocate memory on the stack
# instead of on the heap to minimize memory fragmentation
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90 2008-08-18 00:39:47 UTC (rev 12670)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90 2008-08-18 02:29:55 UTC (rev 12671)
@@ -1931,7 +1931,7 @@
Qmu=600.0d0
Qkappa = 57827.0d0
else
- call exit_mpi_without_rank('Invalid idoubling flag in attenuation_model_1D_prem from get_model()')
+ call exit_mpi_without_rank('Invalid idoubling flag in attenuation_model_1D_prem')
endif
end subroutine attenuation_model_1D_PREM
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90 2008-08-18 00:39:47 UTC (rev 12670)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/comp_mass_matrix_one_element.f90 2008-08-18 02:29:55 UTC (rev 12671)
@@ -1,6 +1,8 @@
!! DK DK added this for merged version
+ if(ipass == 2) then
+
! suppress fictitious elements in central cube
! also take into account the fact that array idoubling is not allocated for the outer core
add_contrib_this_element = .true.
@@ -68,3 +70,5 @@
endif ! of exclusion of fictitious inner core elements
+ endif ! of ipass == 2
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90 2008-08-18 00:39:47 UTC (rev 12670)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90 2008-08-18 02:29:55 UTC (rev 12671)
@@ -45,7 +45,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ipass)
implicit none
@@ -288,7 +288,7 @@
! attenuation_simplex_variables
! correct number of spectral elements in each block depending on chunk type
- integer ispec,nspec,nspec_stacey
+ integer ispec,nspec,nspec_stacey,ipass
integer REFERENCE_1D_MODEL,THREE_D_MODEL
@@ -433,7 +433,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ipass)
! add topography without the crustal model
!! DK DK added this for merged version because array idoubling is not allocated in outer core
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-08-18 00:39:47 UTC (rev 12670)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-08-18 02:29:55 UTC (rev 12671)
@@ -567,6 +567,9 @@
character(len=80) kerstr
character(len=40) varstr(maxker)
+! to perform two passes of the whole routine to be able to save memory
+ integer :: ipass
+
! the height at which the central cube is cut
integer :: nz_inf_limit
@@ -586,6 +589,9 @@
real(kind=CUSTOM_REAL) :: normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
real(kind=CUSTOM_REAL) :: normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
+! perform two passes of the whole routine to be able to save memory
+ do ipass = 1,2
+
! attenuation
if(ATTENUATION .and. ATTENUATION_3D) then
T_c_source = AM_V%QT_c_source
@@ -710,7 +716,7 @@
ystore(:,:,:,:) = 0.d0
zstore(:,:,:,:) = 0.d0
- ibool(:,:,:,:) = 0
+ if(ipass == 1) ibool(:,:,:,:) = 0
! initialize boundary arrays
iboun(:,:) = .false.
@@ -719,7 +725,7 @@
!! DK DK added this for merged version
! creating mass matrix in this slice (will be fully assembled in the solver)
- rmass(:) = 0._CUSTOM_REAL
+ if(ipass == 2) rmass(:) = 0._CUSTOM_REAL
if (.not. PATCH_FOR_GORDON_BELL .and. (CASE_3D .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. SUPPRESS_CRUSTAL_MESH)) then
allocate(stretch_tab(2,ner(1)),STAT=ier )
@@ -853,7 +859,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ipass)
!! DK DK added this for merged version
include "comp_mass_matrix_one_element.f90"
@@ -1037,7 +1043,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ipass)
!! DK DK added this for merged version
include "comp_mass_matrix_one_element.f90"
@@ -1184,7 +1190,6 @@
idoubling(ispec) = IFLAG_IN_FICTITIOUS_CUBE
endif
-
! compute several rheological and geometrical properties for this spectral element
call compute_element_properties(ispec,iregion_code,idoubling, &
xstore,ystore,zstore,nspec, &
@@ -1205,7 +1210,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ipass)
!! DK DK added this for merged version
include "comp_mass_matrix_one_element.f90"
@@ -1222,6 +1227,9 @@
! check total number of spectral elements created
if(ispec /= nspec) call exit_MPI(myrank,'ispec should equal nspec')
+! only create global addressing and the MPI buffers in the first pass
+ if(ipass == 1) then
+
locval = 0
ifseg = .false.
xp = 0.d0
@@ -1249,8 +1257,8 @@
! check that number of points found equals theoretical value
if(nglob /= nglob_theor) then
- write(errmsg,*) 'incorrect total number of points found: myrank,nglob,nglob_theor,iregion_code = ',&
- myrank,nglob,nglob_theor,iregion_code
+ write(errmsg,*) 'incorrect total number of points found: myrank,nglob,nglob_theor,ipass,iregion_code = ',&
+ myrank,nglob,nglob_theor,ipass,iregion_code
call exit_MPI(myrank,errmsg)
endif
@@ -1309,29 +1317,10 @@
iregion_code)
#endif
-! create AVS or DX mesh data for the slices
- if(SAVE_MESH_FILES) then
- call create_name_database(prname,myrank,iregion_code)
+! only create mass matrix and save all the final arrays in the second pass
+ else if(ipass == 2) then
- call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
-
- call write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
-
- call write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot, &
-!! DK DK changed for now because array rhostore is not available in v4.1 anymore
-!! DK DK rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
- nspl,rspl,espl,espl2, &
- ELLIPTICITY,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST,REFERENCE_1D_MODEL, &
- RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
- RMIDDLE_CRUST,ROCEAN,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,iregion_code)
-
- call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
- endif
-
-! copy the theoretical number of points
+! copy the theoretical number of points for the second pass
nglob = nglob_theor
! count number of anisotropic elements in current region
@@ -1445,6 +1434,28 @@
!! DK DK shared by the mesher and the solver subroutines at some point
if(ATTENUATION_VAL) call attenuation_save_arrays(iregion_code, AM_V)
+! create AVS or DX mesh data for the slices
+ if(SAVE_MESH_FILES) then
+ call create_name_database(prname,myrank,iregion_code)
+
+ call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
+
+ call write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
+
+ call write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun,ibool, &
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot, &
+!! DK DK changed for now because array rhostore is not available in v4.1 anymore
+!! DK DK rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST,REFERENCE_1D_MODEL, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,iregion_code)
+
+ call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
+ endif
+
! compute volume, bottom and top area of that part of the slice
volume_local = ZERO
area_local_bottom = ZERO
@@ -1501,5 +1512,11 @@
enddo
enddo
+ else
+ stop 'there cannot be more than two passes in mesh creation'
+ endif ! end of test if first or second pass
+
+ enddo ! of loop on ipass = 1,2
+
end subroutine create_regions_mesh
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90 2008-08-18 00:39:47 UTC (rev 12670)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90 2008-08-18 02:29:55 UTC (rev 12671)
@@ -44,7 +44,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ipass)
implicit none
@@ -321,7 +321,7 @@
real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,nspec_stacey),rho_vs(NGLLX,NGLLY,NGLLZ,nspec_stacey)
- integer nspec_ani
+ integer nspec_ani,ipass
! the 21 coefficients for an anisotropic medium in reduced notation
double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
@@ -389,23 +389,30 @@
character(len=80) kerstr
character(len=40) varstr(maxker)
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
xmesh = ZERO
ymesh = ZERO
zmesh = ZERO
- do ia=1,NGNOD
+
+ do ia = 1,NGNOD
xmesh = xmesh + shape3D(ia,i,j,k)*xelm(ia)
ymesh = ymesh + shape3D(ia,i,j,k)*yelm(ia)
zmesh = zmesh + shape3D(ia,i,j,k)*zelm(ia)
enddo
- r = dsqrt(xmesh*xmesh + ymesh*ymesh + zmesh*zmesh)
xstore(i,j,k) = xmesh
ystore(i,j,k) = ymesh
zstore(i,j,k) = zmesh
+!! DK DK only assign material properties in the second pass, because it involves
+!! DK DK a costly interpolation process at grid points in the case of 3D models
+ if(ipass == 2) then
+
+ r = dsqrt(xmesh*xmesh + ymesh*ymesh + zmesh*zmesh)
+
! make sure we are within the right shell in PREM to honor discontinuities
! use small geometrical tolerance
r_prem = r
@@ -998,6 +1005,8 @@
Qmu_store(i,j,k,ispec) = Qmu
endif
+ endif ! of if(ipass == 2)
+
enddo
enddo
enddo
More information about the cig-commits
mailing list