[cig-commits] [commit] devel, master: better format for the output of Roland_Sylvain gravity integrals, added output of the norm of the g vector, and added an exit message if the solver is run with ROLAND_SYLVAIN = .true. (only the mesher is needed in that case, since we compute some integrals on the mesh and then we are done, no need for a time-domain solver) (c33bf42)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:14:17 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit c33bf4291b15482505f1a8d569d6891a531c820c
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed May 7 02:04:35 2014 +0200
better format for the output of Roland_Sylvain gravity integrals, added output of the norm of the g vector, and added an exit message if the solver is run with ROLAND_SYLVAIN = .true. (only the mesher is needed in that case, since we compute some integrals on the mesh and then we are done, no need for a time-domain solver)
>---------------------------------------------------------------
c33bf4291b15482505f1a8d569d6891a531c820c
setup/constants.h.in | 3 +-
src/meshfem3D/compute_coordinates_grid.f90 | 29 ++++---
src/meshfem3D/finalize_mesher.f90 | 121 +++++++++++++++++++++++++----
src/meshfem3D/meshfem3D_par.f90 | 2 +
src/specfem3D/initialize_simulation.f90 | 5 +-
5 files changed, 131 insertions(+), 29 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 49eacd8..164fe1b 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -420,7 +420,8 @@
! 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
+!! DK DK 2 is a fictitious value used to save memory when the ROLAND_SYLVAIN option is off
+ integer, parameter :: NX_OBSERVATION = 2 ! 500
integer, parameter :: NY_OBSERVATION = NX_OBSERVATION
! the code will display sample output values at this particular point as a check
diff --git a/src/meshfem3D/compute_coordinates_grid.f90 b/src/meshfem3D/compute_coordinates_grid.f90
index 7ddd5ab..35626d2 100644
--- a/src/meshfem3D/compute_coordinates_grid.f90
+++ b/src/meshfem3D/compute_coordinates_grid.f90
@@ -331,7 +331,8 @@
use constants
- use meshfem3D_par, only: myrank,x_observation,y_observation,z_observation,ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES
+ use meshfem3D_par, only: myrank,x_observation,y_observation,z_observation,lon_observation,lat_observation, &
+ ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES
use meshfem3D_models_par, only: nspl,rspl,espl,espl2,ibathy_topo
@@ -420,20 +421,26 @@
! 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)
+ ! 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
+ ! compute elevation at current point
+ if(TOPOGRAPHY) then
call get_topo_bathy(lat,lon,elevation,ibathy_topo)
+ else
+ elevation = ZERO
+ endif
- ! 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
+ ! 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
+ lon_observation(ix,iy,ichunk) = lon
+ lat_observation(ix,iy,ichunk) = lat
+ endif
+ ! add topography
+ if(TOPOGRAPHY) then
! non-dimensionalize the elevation, which is in meters
elevation = elevation / R_EARTH
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index e04623b..5870b77 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -128,28 +128,28 @@
write(IMAIN,*) 'maxval of norm of g vector on whole observation surface = ',maxval(temporary_array_for_sum),' m.s-2'
write(IMAIN,*)
- write(IMAIN,*) 'minval of abs(G_xx) on whole observation surface = ',minval(abs(G_xx)),' Eotvos'
- write(IMAIN,*) 'maxval of abs(G_xx) on whole observation surface = ',maxval(abs(G_xx)),' Eotvos'
+ write(IMAIN,*) 'minval of G_xx on whole observation surface = ',minval(G_xx),' Eotvos'
+ write(IMAIN,*) 'maxval of G_xx on whole observation surface = ',maxval(G_xx),' Eotvos'
write(IMAIN,*)
- write(IMAIN,*) 'minval of abs(G_yy) on whole observation surface = ',minval(abs(G_yy)),' Eotvos'
- write(IMAIN,*) 'maxval of abs(G_yy) on whole observation surface = ',maxval(abs(G_yy)),' Eotvos'
+ write(IMAIN,*) 'minval of G_yy on whole observation surface = ',minval(G_yy),' Eotvos'
+ write(IMAIN,*) 'maxval of G_yy on whole observation surface = ',maxval(G_yy),' Eotvos'
write(IMAIN,*)
- write(IMAIN,*) 'minval of abs(G_zz) on whole observation surface = ',minval(abs(G_zz)),' Eotvos'
- write(IMAIN,*) 'maxval of abs(G_zz) on whole observation surface = ',maxval(abs(G_zz)),' Eotvos'
+ write(IMAIN,*) 'minval of G_zz on whole observation surface = ',minval(G_zz),' Eotvos'
+ write(IMAIN,*) 'maxval of G_zz on whole observation surface = ',maxval(G_zz),' Eotvos'
write(IMAIN,*)
- write(IMAIN,*) 'minval of abs(G_xy) on whole observation surface = ',minval(abs(G_xy)),' Eotvos'
- write(IMAIN,*) 'maxval of abs(G_xy) on whole observation surface = ',maxval(abs(G_xy)),' Eotvos'
+ write(IMAIN,*) 'minval of G_xy on whole observation surface = ',minval(G_xy),' Eotvos'
+ write(IMAIN,*) 'maxval of G_xy on whole observation surface = ',maxval(G_xy),' Eotvos'
write(IMAIN,*)
- write(IMAIN,*) 'minval of abs(G_xz) on whole observation surface = ',minval(abs(G_xz)),' Eotvos'
- write(IMAIN,*) 'maxval of abs(G_xz) on whole observation surface = ',maxval(abs(G_xz)),' Eotvos'
+ write(IMAIN,*) 'minval of G_xz on whole observation surface = ',minval(G_xz),' Eotvos'
+ write(IMAIN,*) 'maxval of G_xz on whole observation surface = ',maxval(G_xz),' Eotvos'
write(IMAIN,*)
- write(IMAIN,*) 'minval of abs(G_yz) on whole observation surface = ',minval(abs(G_yz)),' Eotvos'
- write(IMAIN,*) 'maxval of abs(G_yz) on whole observation surface = ',maxval(abs(G_yz)),' Eotvos'
+ write(IMAIN,*) 'minval of G_yz on whole observation surface = ',minval(G_yz),' Eotvos'
+ write(IMAIN,*) 'maxval of G_yz on whole observation surface = ',maxval(G_yz),' Eotvos'
write(IMAIN,*)
write(IMAIN,*) 'Minval and maxval of trace of G, which in principle should be zero:'
@@ -196,14 +196,103 @@
write(IMAIN,*) 'computed G_yz = ',G_yz(ixr,iyr,ichunkr),' Eotvos'
! for future GMT display
- open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_g_and_G_for_GMT.txt',status='unknown',action='write')
! loop on all the chunks and then on all the observation nodes in each chunk
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_g_x_for_GMT.txt',status='unknown',action='write')
do ichunkval = 1,NCHUNKS_MAX
do iyval = 1,NY_OBSERVATION
do ixval = 1,NX_OBSERVATION
- write(IOUT,*) g_x(ixval,iyval,ichunkval),g_y(ixval,iyval,ichunkval),g_z(ixval,iyval,ichunkval), &
- G_xx(ixval,iyval,ichunkval),G_yy(ixval,iyval,ichunkval),G_zz(ixval,iyval,ichunkval), &
- G_xy(ixval,iyval,ichunkval),G_xz(ixval,iyval,ichunkval),G_yz(ixval,iyval,ichunkval)
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),g_x(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_g_y_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),g_y(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_g_z_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),g_z(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_norm_of_g_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval), &
+ sqrt(g_x(ixval,iyval,ichunkval)**2 + g_y(ixval,iyval,ichunkval)**2 + g_z(ixval,iyval,ichunkval)**2)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_G_xx_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),G_xx(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_G_yy_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),G_yy(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_G_zz_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),G_zz(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_G_xy_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),G_xy(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_G_xz_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),G_xz(ixval,iyval,ichunkval)
+ enddo
+ enddo
+ enddo
+ close(unit=IOUT)
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_G_yz_for_GMT.txt',status='unknown',action='write')
+ do ichunkval = 1,NCHUNKS_MAX
+ do iyval = 1,NY_OBSERVATION
+ do ixval = 1,NX_OBSERVATION
+ write(IOUT,*) lon_observation(ixval,iyval,ichunkval),lat_observation(ixval,iyval,ichunkval),G_yz(ixval,iyval,ichunkval)
enddo
enddo
enddo
diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90
index dcdc493..c3dfc76 100644
--- a/src/meshfem3D/meshfem3D_par.f90
+++ b/src/meshfem3D/meshfem3D_par.f90
@@ -167,6 +167,8 @@
equivalence(y_observation,y_observation1D)
equivalence(z_observation,z_observation1D)
+ double precision, dimension(NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) :: lon_observation,lat_observation
+
! arrays containing the computed fields for Roland_Sylvain integrals at the observation points
! the 1D equivalenced versions are for the FORCE_VECTORIZATION version of the loops
double precision, dimension(NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) :: g_x,g_y,g_z,G_xx,G_yy,G_zz,G_xy,G_xz,G_yz, &
diff --git a/src/specfem3D/initialize_simulation.f90 b/src/specfem3D/initialize_simulation.f90
index fa62b7c..e9a6518 100644
--- a/src/specfem3D/initialize_simulation.f90
+++ b/src/specfem3D/initialize_simulation.f90
@@ -44,6 +44,9 @@
call world_size(sizeprocs)
call world_rank(myrank)
+!! DK DK for Roland_Sylvain
+ if(ROLAND_SYLVAIN) call exit_MPI(myrank,'no need to run the solver to compute Roland_Sylvain integrals, only the mesher')
+
if (myrank == 0) then
! read the parameter file and compute additional parameters
call read_compute_parameters()
@@ -55,7 +58,7 @@
! check that the code is running with the requested nb of processes
if( sizeprocs /= NPROCTOT ) then
print*,'error: rank ',myrank,' - wrong number of MPI processes',sizeprocs,NPROCTOT
- call exit_MPI(myrank,'wrong number of MPI processes(initialization specfem)')
+ call exit_MPI(myrank,'wrong number of MPI processes in the initialization of SPECFEM')
endif
! synchronizes processes
More information about the CIG-COMMITS
mailing list