Skip to content

Commit

Permalink
Merge pull request #72 from komatits/devel
Browse files Browse the repository at this point in the history 
switched from a single observation point to a whole surface around the E...
  • Loading branch information
@komatits
komatits committed May 5, 2014
2 parents a2b53e8 + 47332f7 commit 7be81d6
Show file tree
Hide file tree
Showing 13 changed files with 542 additions and 202 deletions.
49 changes: 28 additions & 21 deletions setup/constants.h.in
View file
Expand Up @@ -91,6 +91,8 @@
! average density in the full Earth to normalize equation
double precision, parameter :: RHOAV = 5514.3d0

! maximum number of chunks (full sphere)
integer, parameter :: NCHUNKS_MAX = 6

!!-----------------------------------------------------------
!!
Expand Down Expand Up @@ -402,30 +404,35 @@
!! DK DK for Roland_Sylvain
logical, parameter :: ROLAND_SYLVAIN = .false.

! altitude of the observation point
double precision, parameter :: altitude_of_observation_point = 255.d0 ! in km
! altitude of the observation points in km
double precision, parameter :: altitude_of_observation_points = 255.d0

! standard gravity at the surface of the Earth
double precision, parameter :: STANDARD_GRAVITY_EARTH = 9.80665d0 ! in m.s-2
! elevation ratio of the points at which we observe
double precision, parameter :: observation_elevation_ratio = (R_EARTH_KM + altitude_of_observation_points) / R_EARTH_KM

! compute the contribution of the crust only (otherwise by default compute the contribution of the whole Earth)
logical, parameter :: COMPUTE_CRUST_CONTRIB_ONLY = .false.

double precision, parameter :: SQUARE_ROOT_OF_THREE = 1.73205080756887729352d0
! check for negative Jacobians in the calculation of integrals or not
! (can safely be done once to check that the mesh is OK at a given resolution, and then permanently
! turned off in future runs because the mesh does not change)
logical, parameter :: CHECK_FOR_NEGATIVE_JACOBIANS = .true.

! position of the observation point in non-dimensionalized value
! number of points in each horizontal direction of the observation grid of each cubed-sphere chunk
! at the altitude of the observation point
integer, parameter :: NX_OBSERVATION = 500
integer, parameter :: NY_OBSERVATION = NX_OBSERVATION

!! DK DK along X only
! double precision, parameter :: x_observation = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM
! double precision, parameter :: y_observation = 0.d0
! double precision, parameter :: z_observation = 0.d0
! the code will display sample output values at this particular point as a check
integer, parameter :: ixr = NX_OBSERVATION / 3
integer, parameter :: iyr = NY_OBSERVATION / 4
integer, parameter :: ichunkr = 3

!! DK DK along diagonal
! double precision, parameter :: x_observation = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE
! double precision, parameter :: y_observation = x_observation
! double precision, parameter :: z_observation = x_observation
! how often (every how many spectral elements computed) we print a timestamp to monitor the behavior of the code
integer, parameter :: NSPEC_DISPLAY_INTERVAL = 100

!! DK DK point randomly chosen in space, not at the right elevation we want
double precision, parameter :: x_observation = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE
double precision, parameter :: y_observation = x_observation * 1.12d0
double precision, parameter :: z_observation = x_observation * 1.38d0
! for the FORCE_VECTORIZATION 1D version of some loops
integer, parameter :: NTOTAL_OBSERVATION = NX_OBSERVATION * NY_OBSERVATION * NCHUNKS_MAX

!------------------------------------------------------
!----------- do not modify anything below -------------
Expand Down Expand Up @@ -471,6 +478,9 @@
! double precision, parameter :: GRAV = 6.6723d-11
double precision, parameter :: GRAV = 6.67384d-11

! standard gravity at the surface of the Earth
double precision, parameter :: STANDARD_GRAVITY_EARTH = 9.80665d0 ! in m.s-2

! a few useful constants
double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
double precision, parameter :: ONE_HALF = HALF, ONE_FOURTH = 0.25d0, ONE_EIGHTH = 0.125d0, ONE_SIXTEENTH = 0.0625d0
Expand Down Expand Up @@ -505,9 +515,6 @@
! Mountains can create signals of several hundred Eotvos.
double precision, parameter :: SI_UNITS_TO_EOTVOS = 1.d+9

! maximum number of chunks (full sphere)
integer, parameter :: NCHUNKS_MAX = 6

! define block type based upon chunk number (between 1 and 6)
! do not change this numbering, chunk AB must be number 1 for central cube
integer, parameter :: CHUNK_AB = 1
Expand Down
161 changes: 148 additions & 13 deletions src/meshfem3D/compute_coordinates_grid.f90
View file
Expand Up @@ -57,8 +57,7 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, &

double precision, dimension(NDIM) :: vector_ori,vector_rotated

double precision :: ratio_xi, ratio_eta, fact_xi, fact_eta, &
fact_xi_,fact_eta_
double precision :: ratio_xi, ratio_eta, fact_xi, fact_eta, fact_xi_,fact_eta_

! this to avoid compilation warnings
x_=0
Expand All @@ -82,8 +81,8 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, &
! uncomment the corresponding lines in the else condition of this if statement too.
! note that the ratio bigger_edge_size/smaller_edge_size for the surface mesh is a bit higher (1.43 vs 1.41)

! fact_xi_= (3.d0*fact_xi+4.d0*fact_xi_)/7.d0
! fact_eta_= (3.d0*fact_eta+4.d0*fact_eta_)/7.d0
! fact_xi_= (3.d0*fact_xi+4.d0*fact_xi_)/7.d0
! fact_eta_= (3.d0*fact_eta+4.d0*fact_eta_)/7.d0

xi = PI_OVER_TWO*fact_xi
eta = PI_OVER_TWO*fact_eta
Expand Down Expand Up @@ -134,15 +133,15 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, &
case default
stop 'incorrect chunk number in compute_coord_main_mesh'
end select
! write(IMAIN,*) x,' ',y,' ',z

else

! uncomment the following lines to have more regular surface mesh (better aspect ratio for each element)
! note that the ratio bigger_edge_size/smaller_edge_size for the surface mesh is a bit higher (1.43 vs 1.41)
! ratio_xi = ((iproc_xi + offset_x(ignod)/dble(NEX_PER_PROC_XI))/dble(NPROC_XI))*tan(ANGULAR_WIDTH_XI_RAD/2.d0)
! x_ = 2.d0*ratio_xi-tan(ANGULAR_WIDTH_XI_RAD/2.d0)
! ratio_eta = ((iproc_eta + offset_y(ignod)/dble(NEX_PER_PROC_ETA))/dble(NPROC_ETA))*tan(ANGULAR_WIDTH_ETA_RAD/2.d0)
! y_ = 2.d0*ratio_eta-tan(ANGULAR_WIDTH_ETA_RAD/2.d0)
! ratio_xi = ((iproc_xi + offset_x(ignod)/dble(NEX_PER_PROC_XI))/dble(NPROC_XI))*tan(ANGULAR_WIDTH_XI_RAD/2.d0)
! x_ = 2.d0*ratio_xi-tan(ANGULAR_WIDTH_XI_RAD/2.d0)
! ratio_eta = ((iproc_eta + offset_y(ignod)/dble(NEX_PER_PROC_ETA))/dble(NPROC_ETA))*tan(ANGULAR_WIDTH_ETA_RAD/2.d0)
! y_ = 2.d0*ratio_eta-tan(ANGULAR_WIDTH_ETA_RAD/2.d0)

ratio_xi = ((iproc_xi + offset_x(ignod)/dble(NEX_PER_PROC_XI))/dble(NPROC_XI))
x = 2.d0*ratio_xi-1
Expand All @@ -155,8 +154,8 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, &

! uncomment the following lines to have more regular surface mesh (better aspect ratio for each element)
! note that the ratio bigger_edge_size/smaller_edge_size for the surface mesh is a bit higher (1.43 vs 1.41)
! x= (3.d0*x_+4.d0*x)/7.d0
! y= (3.d0*y_+4.d0*y)/7.d0
! x= (3.d0*x_+4.d0*x)/7.d0
! y= (3.d0*y_+4.d0*y)/7.d0

gamma = ONE / sqrt(ONE + x*x + y*y)

Expand Down Expand Up @@ -272,12 +271,12 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, &

endif
enddo
! if(ilayer == NUMBER_OF_MESH_LAYERS .and. INCLUDE_CENTRAL_CUBE) write(IMAIN,*)

end subroutine compute_coord_main_mesh

!---------------------------------------------------------------------------

!! DK DK create value of arrays xgrid ygrid and zgrid in the central cube without storing them
! create value of arrays xgrid ygrid and zgrid in the central cube without storing them

subroutine compute_coord_central_cube(ix,iy,iz, &
xgrid_central_cube,ygrid_central_cube,zgrid_central_cube, &
Expand Down Expand Up @@ -321,3 +320,139 @@ subroutine compute_coord_central_cube(ix,iy,iz, &

end subroutine compute_coord_central_cube


!---------------------------------------------------------------------------

! create observation grid (surface) on which we compute the Roland_Sylvain integrals.
! each processor has a copy of the whole grid and sums its contribution there,
! and at the end of the process an MPI reduction is performed to compute the global sum

subroutine compute_observation_surface()

use constants

use meshfem3D_par, only: myrank,x_observation,y_observation,z_observation,ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES

use meshfem3D_models_par, only: nspl,rspl,espl,espl2,ibathy_topo

implicit none

! local variables
integer :: ix,iy,ichunk

double precision :: gamma,x,y,rgt
double precision :: x_top,y_top,z_top
double precision :: ratio_xi, ratio_eta

double precision :: r,lat,lon,elevation

! the observation surface is always above the whole Earth, thus each mesh chunk has a size of PI / 2 in each direction
double precision, parameter :: ANGULAR_WIDTH_XI_RAD = PI_OVER_TWO, ANGULAR_WIDTH_ETA_RAD = PI_OVER_TWO

! for future GMT display
if(myrank == 0) open(unit=IOUT,file=trim(OUTPUT_FILES)//'/observation_grid_long_lat_topo_for_GMT.txt', &
status='unknown',action='write')

! loop on all the chunks and then on all the observation nodes in each chunk
do ichunk = 1,NCHUNKS_MAX
do iy = 1,NY_OBSERVATION
do ix = 1,NX_OBSERVATION

ratio_xi = dble(ix - 1) / dble(NX_OBSERVATION - 1)
x = 2.d0*ratio_xi-1

ratio_eta = dble(iy - 1) / dble(NY_OBSERVATION - 1)
y = 2.d0*ratio_eta-1

x = tan((ANGULAR_WIDTH_XI_RAD/2.d0) * x)
y = tan((ANGULAR_WIDTH_ETA_RAD/2.d0) * y)

gamma = ONE / sqrt(ONE + x*x + y*y)

! first compute the position of the points exactly at the free surface of the Earth (without the oceans)
! keeping in mind that the code non-dimensionalizes the radius of the spherical Earth to one
rgt = R_UNIT_SPHERE*gamma

! define the mesh points on the top and the bottom in the six regions of the cubed shpere
select case (ichunk)

case(CHUNK_AB)

x_top = -y*rgt
y_top = x*rgt
z_top = rgt

case(CHUNK_AB_ANTIPODE)

x_top = -y*rgt
y_top = -x*rgt
z_top = -rgt

case(CHUNK_AC)

x_top = -y*rgt
y_top = -rgt
z_top = x*rgt

case(CHUNK_AC_ANTIPODE)

x_top = -y*rgt
y_top = rgt
z_top = -x*rgt

case(CHUNK_BC)

x_top = -rgt
y_top = y*rgt
z_top = x*rgt

case(CHUNK_BC_ANTIPODE)

x_top = rgt
y_top = -y*rgt
z_top = x*rgt

case default
stop 'incorrect chunk number in compute_coord_main_mesh'

end select

! add ellipticity
if(ELLIPTICITY) call get_ellipticity_single_point(x_top,y_top,z_top,nspl,rspl,espl,espl2)

! add topography
if(TOPOGRAPHY) then
! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees)
call xyz_2_rlatlon_dble(x_top,y_top,z_top,r,lat,lon)

! compute elevation at current point
call get_topo_bathy(lat,lon,elevation,ibathy_topo)

! store the values obtained for future display with GMT
if(myrank == 0) then
if( lon > 180.0d0 ) lon = lon - 360.0d0
write(IOUT,*) lon,lat,elevation
endif

! non-dimensionalize the elevation, which is in meters
elevation = elevation / R_EARTH

x_top = x_top*(ONE + elevation)
y_top = y_top*(ONE + elevation)
z_top = z_top*(ONE + elevation)
endif

! compute the position of the point
x_observation(ix,iy,ichunk) = x_top * observation_elevation_ratio
y_observation(ix,iy,ichunk) = y_top * observation_elevation_ratio
z_observation(ix,iy,ichunk) = z_top * observation_elevation_ratio

enddo
enddo
enddo

! for future GMT display
if(myrank == 0) close(unit=IOUT)

end subroutine compute_observation_surface

0 comments on commit 7be81d6

Please sign in to comment.