[cig-commits] r12605 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: DATA src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Aug 11 13:50:58 PDT 2008
Author: dkomati1
Date: 2008-08-11 13:50:57 -0700 (Mon, 11 Aug 2008)
New Revision: 12605
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem1.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem2.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem1.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem2.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_main.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/memory_eval.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90
Log:
Put the implementation of the oceans back in the merged version.
Fixed a small issue (incorrect flag) when computing the local height of oceans in the case of a 3D Earth.
Added stop statements if GRAVITY, ROTATION or SAVE_FORWARD are true or if SIMULATION_TYPE > 1 in DATA/Par_file because for now the merged version can only handle a high-frequency forward problem.
Removed some useless white spaces.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file 2008-08-11 20:50:57 UTC (rev 12605)
@@ -15,7 +15,7 @@
# number of elements at the surface along the two sides of the first chunk
# (must be multiple of 16 and also of 8 * multiple of NPROC below)
-NEX_XI = 128 # 144 288 432 576 720 864 1008 1152 1296 1440
+NEX_XI = 128
NEX_ETA = 128
# number of MPI processors along the two sides of the first chunk
@@ -31,12 +31,12 @@
# fully 3D models:
# transversely_isotropic_prem_plus_3D_crust_2.0, 3D_anisotropic, 3D_attenuation,
# s20rts, s362ani, s362iso, s362wmani, s362ani_prem, s29ea, s29ea,sea99_jp3d1994,sea99,jp3d1994
-MODEL = 1D_transversely_isotropic_prem_onecrust
+MODEL = 1D_transversely_isotropic_prem
# parameters describing the Earth model
-OCEANS = .false.
-ELLIPTICITY = .false.
-TOPOGRAPHY = .false.
+OCEANS = .true.
+ELLIPTICITY = .true.
+TOPOGRAPHY = .true.
GRAVITY = .false.
ROTATION = .false.
ATTENUATION = .false.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem1.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem1.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem1.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -19,5 +19,6 @@
iboolcorner_crust_mantle,iboolcorner_outer_core,iboolcorner_inner_core, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
npoin2D_faces_outer_core,npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
- npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core)
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core,rmass_ocean_load, &
+ normal_top_crust_mantle,ibelm_top_crust_mantle)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem2.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem2.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_meshfem2.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -19,5 +19,6 @@
iboolcorner_crust_mantle,iboolcorner_outer_core,iboolcorner_inner_core, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
npoin2D_faces_outer_core,npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
- npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core)
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core,rmass_ocean_load, &
+ normal_top_crust_mantle,ibelm_top_crust_mantle)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem1.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem1.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem1.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -24,5 +24,6 @@
iboolcorner_crust_mantle,iboolcorner_outer_core,iboolcorner_inner_core, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
npoin2D_faces_outer_core,npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
- npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core)
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core,rmass_ocean_load, &
+ normal_top_crust_mantle,ibelm_top_crust_mantle)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem2.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem2.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/call_specfem2.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -24,5 +24,6 @@
iboolcorner_crust_mantle,iboolcorner_outer_core,iboolcorner_inner_core, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
npoin2D_faces_outer_core,npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
- npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core)
+ npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core,rmass_ocean_load, &
+ normal_top_crust_mantle,ibelm_top_crust_mantle)
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-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -1364,8 +1364,8 @@
iglobnum=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
-! compute local height of oceans
- if(ISOTROPIC_3D_MANTLE) then
+! if 3D Earth, compute local height of oceans
+ if(CASE_3D) then
! get coordinates of current point
xval = xstore(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
@@ -1396,6 +1396,7 @@
endif
else
+! if 1D Earth, use oceans of constant thickness everywhere
height_oceans = THICKNESS_OCEANS_PREM
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_main.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_main.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_main.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -15,6 +15,10 @@
npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
npoin2D_xi_inner_core,npoin2D_eta_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_CM) :: normal_top_crust_mantle
+ integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
+
! number of elements on the boundaries
integer :: nspec2D_xmin_inner_core,nspec2D_xmax_inner_core,nspec2D_ymin_inner_core,nspec2D_ymax_inner_core
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -36,6 +36,8 @@
double precision, dimension(NSPEC_CRUST_MANTLE * NGLLX * NGLLY * NGLLZ) :: xp,yp,zp
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_CM) :: normal_top_crust_mantle
+ integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
integer, dimension(NGLOB_CRUST_MANTLE) :: mask_ibool
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: copy_ibool_ori
@@ -99,7 +101,6 @@
integer, dimension(NSPEC2DMAX_XMIN_XMAX_CM) :: ibelm_xmin_crust_mantle,ibelm_xmax_crust_mantle
integer, dimension(NSPEC2DMAX_YMIN_YMAX_CM) :: ibelm_ymin_crust_mantle,ibelm_ymax_crust_mantle
integer, dimension(NSPEC2D_BOTTOM_CM) :: ibelm_bottom_crust_mantle
- integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX_CM) :: &
jacobian2D_xmin_crust_mantle,jacobian2D_xmax_crust_mantle
@@ -111,7 +112,6 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX_CM) :: normal_xmin_crust_mantle,normal_xmax_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2DMAX_YMIN_YMAX_CM) :: normal_ymin_crust_mantle,normal_ymax_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM_CM) :: normal_bottom_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_CM) :: normal_top_crust_mantle
integer, dimension(NSPEC2DMAX_XMIN_XMAX_OC) :: ibelm_xmin_outer_core,ibelm_xmax_outer_core
integer, dimension(NSPEC2DMAX_YMIN_YMAX_OC) :: ibelm_ymin_outer_core,ibelm_ymax_outer_core
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-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -868,7 +868,7 @@
! store this as well for anisotropic elements in the mantle
if(iregion_code == IREGION_CRUST_MANTLE .and. NSPECMAX_TISO_MANTLE > 1) then
-! only store the anisotropic layers, because these arrays have purposely been
+! only store the anisotropic layers, because these arrays have purposely been
! declared with a reduced size to save memory
if(ispec <= NSPECMAX_TISO_MANTLE) then
kappahstore(i,j,k,ispec) = sngl(rho*(vph*vph - 4.d0*vsh*vsh/3.d0))
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -207,7 +207,7 @@
!! DK DK Section 5.6 about this
module dyn_array
!---------------------------------------------------------------------
-! Module containing definitions needed to dynamically allocate the values of an array
+! Module containing definitions needed to dynamically allocate the values of an array
!---------------------------------------------------------------------
include "constants.h"
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/memory_eval.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/memory_eval.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -299,6 +299,9 @@
! rmass_ocean_load
static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
+! updated_dof_ocean_load
+ static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*dble(SIZE_LOGICAL)
+
endif
! add arrays used to save strain for attenuation or for adjoint runs
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -152,6 +152,10 @@
call read_value_logical(SAVE_FORWARD, 'solver.SAVE_FORWARD')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+ if(SIMULATION_TYPE > 1) stop 'SIMULATION_TYPE > 1 not implemented in the reduced merged version yet'
+
+ if(SAVE_FORWARD) stop 'SAVE_FORWARD not implemented in the reduced merged version yet'
+
call read_value_integer(NCHUNKS, 'mesher.NCHUNKS')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
if(NCHUNKS /= 1 .and. NCHUNKS /= 2 .and. NCHUNKS /= 3 .and. NCHUNKS /= 6) &
@@ -765,7 +769,7 @@
! take a 5% safety margin on the maximum stable time step
! which was obtained by trial and error
- DT = DT * (1.d0 - 0.05d0)
+!!!!!!!!!!!!!!!!!! DT = DT * (1.d0 - 0.05d0)
call read_value_logical(OCEANS, 'model.OCEANS')
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
@@ -789,6 +793,10 @@
if(ATTENUATION_3D .and. .not. ATTENUATION) stop 'need ATTENUATION to use ATTENUATION_3D'
+ if(GRAVITY) stop 'GRAVITY not implemented in the reduced merged version yet because useless at high frequency'
+
+ if(ROTATION) stop 'ROTATION not implemented in the reduced merged version yet because useless at high frequency'
+
! radii in PREM or IASP91
! and normalized density at fluid-solid interface on fluid size for coupling
! ROCEAN: radius of the ocean (m)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90 2008-08-11 20:36:00 UTC (rev 12604)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90 2008-08-11 20:50:57 UTC (rev 12605)
@@ -110,7 +110,14 @@
! additional mass matrix for ocean load
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_CM) :: normal_top_crust_mantle
+ integer, dimension(NSPEC2D_TOP_CM) :: ibelm_top_crust_mantle
+! flag to mask ocean-bottom degrees of freedom for ocean load
+ logical, dimension(NGLOB_CRUST_MANTLE_OCEANS) :: updated_dof_ocean_load
+
+ real(kind=CUSTOM_REAL) additional_term,force_normal_comp
+
! arrays to couple with the fluid regions by pointwise matching
integer, dimension(NSPEC2D_BOTTOM_OC) :: ibelm_bottom_outer_core
integer, dimension(NSPEC2D_TOP_OC) :: ibelm_top_outer_core
@@ -783,9 +790,6 @@
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROCTOT) call exit_MPI(myrank,'wrong number of MPI processes')
-!! DK DK added this
- if(OCEANS) call exit_MPI(myrank,'DK DK je crois que j ai enleve les oceans par erreur, les remettre')
-
! check that the code has been compiled with the right values
if (NSPEC_computed(IREGION_CRUST_MANTLE) /= NSPEC_CRUST_MANTLE) then
write(IMAIN,*) NSPEC_computed(IREGION_CRUST_MANTLE),NSPEC_CRUST_MANTLE
@@ -2090,6 +2094,56 @@
accel_crust_mantle(3,i) = accel_crust_mantle(3,i)*rmass_crust_mantle(i)
enddo
+ if(OCEANS) then
+
+! initialize the updates
+ updated_dof_ocean_load(:) = .false.
+
+! for surface elements exactly at the top of the crust (ocean bottom)
+ do ispec2D = 1,NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+
+ ispec = ibelm_top_crust_mantle(ispec2D)
+
+! only for DOFs exactly at the top of the crust (ocean bottom)
+ k = NGLLZ
+
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+! get global point number
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+
+! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
+
+! get normal
+ nx = normal_top_crust_mantle(1,i,j,ispec2D)
+ ny = normal_top_crust_mantle(2,i,j,ispec2D)
+ nz = normal_top_crust_mantle(3,i,j,ispec2D)
+
+! make updated component of right-hand side
+! we divide by rmass_crust_mantle() which is 1 / M
+! we use the total force which includes the Coriolis term above
+ force_normal_comp = (accel_crust_mantle(1,iglob)*nx + &
+ accel_crust_mantle(2,iglob)*ny + &
+ accel_crust_mantle(3,iglob)*nz) / rmass_crust_mantle(iglob)
+
+ additional_term = (rmass_ocean_load(iglob) - rmass_crust_mantle(iglob)) * force_normal_comp
+
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term * nx
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term * ny
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term * nz
+
+! done with this point
+ updated_dof_ocean_load(iglob) = .true.
+
+ endif
+
+ enddo
+ enddo
+ enddo
+ endif
+
do i=1,NGLOB_CRUST_MANTLE
veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
enddo
More information about the cig-commits
mailing list