[cig-commits] [commit] devel, master: fixed physical units of Roland_Sylvain integrals (cb1e8f4)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:13:10 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit cb1e8f4563a385d39999cfd97a12ddc693b8ef77
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu May 1 14:42:44 2014 +0200
fixed physical units of Roland_Sylvain integrals
>---------------------------------------------------------------
cb1e8f4563a385d39999cfd97a12ddc693b8ef77
setup/constants.h.in | 52 +++++++++++++++++++++----
src/meshfem3D/compute_volumes_and_areas.f90 | 29 ++++----------
src/meshfem3D/finalize_mesher.f90 | 59 +++++++++++++++++++----------
3 files changed, 92 insertions(+), 48 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 655f912..79fa1b9 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -48,9 +48,6 @@
!*********************************************************************************************************
-!! DK DK for Roland_Sylvain
- logical, parameter :: ROLAND_SYLVAIN = .false.
-
! if files on a local path on each node are also seen as global with same path
! set to .true. typically on a machine with a common (shared) file system, e.g. LUSTRE, GPFS or NFS-mounted /home
! set to .false. typically on a cluster of nodes with local disks used to store the mesh (not very common these days)
@@ -88,6 +85,9 @@
! uncomment line below for PREM with oceans
! double precision, parameter :: R_EARTH = 6368000.d0
+! radius of the Earth in km
+ double precision, parameter :: R_EARTH_KM = R_EARTH / 1000.d0
+
! average density in the full Earth to normalize equation
double precision, parameter :: RHOAV = 5514.3d0
@@ -393,6 +393,41 @@
integer, parameter :: NSTEP_FOR_BENCHMARK = 300
logical, parameter :: SET_INITIAL_FIELD_TO_1_IN_BENCH = .true.
+!!-----------------------------------------------------------
+!!
+!! for Roland_Sylvain integrals
+!!
+!!-----------------------------------------------------------
+
+!! 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
+
+! standard gravity at the surface of the Earth
+ double precision, parameter :: STANDARD_GRAVITY_EARTH = 9.80665d0 ! in m.s-2
+
+ double precision, parameter :: SQUARE_ROOT_OF_THREE = 1.73205080756887729352d0
+
+! position of the observation point in non-dimensionalized value
+
+!! DK DK along X only
+! double precision, parameter :: xobs = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM
+! double precision, parameter :: yobs = 0.d0
+! double precision, parameter :: zobs = 0.d0
+
+!! DK DK along diagonal
+! double precision, parameter :: xobs = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE
+! double precision, parameter :: yobs = xobs
+! double precision, parameter :: zobs = xobs
+
+!! DK DK point randomly chosen in space, not at the right elevation we want
+ double precision, parameter :: xobs = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE
+ double precision, parameter :: yobs = xobs * 1.12d0
+ double precision, parameter :: zobs = xobs * 1.38d0
+
+
!------------------------------------------------------
!----------- do not modify anything below -------------
!------------------------------------------------------
@@ -430,8 +465,12 @@
integer, parameter :: MIDY = (NGLLY+1)/2
integer, parameter :: MIDZ = (NGLLZ+1)/2
-! gravitational constant in S.I. units i.e. in m3 kg-1 s-2
- double precision, parameter :: GRAV = 6.6723d-11
+! gravitational constant in S.I. units i.e. in m3 kg-1 s-2, or equivalently in N.(m/kg)^2
+!! DK DK April 2014: switched to the 2010 Committee on Data for Science and Technology (CODATA) recommended value
+!! DK DK see e.g. http://www.physics.nist.gov/cgi-bin/cuu/Value?bg
+!! DK DK and http://en.wikipedia.org/wiki/Gravitational_constant
+! double precision, parameter :: GRAV = 6.6723d-11
+ double precision, parameter :: GRAV = 6.67384d-11
! a few useful constants
double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
@@ -454,9 +493,6 @@
! normalized radius of free surface
double precision, parameter :: R_UNIT_SPHERE = ONE
-! radius of the Earth in km
- double precision, parameter :: R_EARTH_KM = R_EARTH / 1000.d0
-
! fixed thickness of 3 km for PREM oceans
double precision, parameter :: THICKNESS_OCEANS_PREM = 3000.d0 / R_EARTH
diff --git a/src/meshfem3D/compute_volumes_and_areas.f90 b/src/meshfem3D/compute_volumes_and_areas.f90
index adcf05b..050408c 100644
--- a/src/meshfem3D/compute_volumes_and_areas.f90
+++ b/src/meshfem3D/compute_volumes_and_areas.f90
@@ -294,27 +294,13 @@
double precision, dimension(9) :: Roland_Sylvain_int_local,Roland_Sylvain_int_total_region,elemental_contribution
! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized
- double precision, parameter :: non_dimensionalizing_factor = RHOAV*R_EARTH**3
-
-!! DK DK position of the observation point in non-dimensional value
- double precision, parameter :: altitude_of_observation_point = 255.d0 ! in km
-
- double precision, parameter :: SQUARE_ROOT_OF_THREE = 1.732050808d0
-
-!! DK DK along X only
-! double precision, parameter :: xobs = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM
-! double precision, parameter :: yobs = 0.d0
-! double precision, parameter :: zobs = 0.d0
-
-!! DK DK along diagonal
- double precision, parameter :: xobs = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE
- double precision, parameter :: yobs = xobs
- double precision, parameter :: zobs = xobs
+ ! for the gravity vector force, a distance is involved in the dimensions
+ double precision, parameter :: nondimensionalizing_factor_gi = RHOAV * R_EARTH
+ ! for the second-order gravity tensor, no distance is involved in the dimensions
+ double precision, parameter :: nondimensionalizing_factor_Gij = RHOAV
-!! DK DK point randomly chosen in space, not at the right elevation we want
-! double precision, parameter :: xobs = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE
-! double precision, parameter :: yobs = xobs * 1.12d0
-! double precision, parameter :: zobs = xobs * 1.38d0
+ double precision, parameter :: scaling_factor_gi = GRAV * nondimensionalizing_factor_gi
+ double precision, parameter :: scaling_factor_Gij = GRAV * nondimensionalizing_factor_Gij
! initializes
Roland_Sylvain_int_local(:) = ZERO
@@ -399,7 +385,8 @@
! multiply by the gravitational constant in S.I. units i.e. in m3 kg-1 s-2
! and also take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized
- Roland_Sylvain_int_local(:) = Roland_Sylvain_int_local(:) * GRAV * non_dimensionalizing_factor
+ Roland_Sylvain_int_local(1:3) = Roland_Sylvain_int_local(1:3) * scaling_factor_gi
+ Roland_Sylvain_int_local(4:9) = Roland_Sylvain_int_local(4:9) * scaling_factor_Gij
! use an MPI reduction to compute the total value of the integral
Roland_Sylvain_int_total_region(:) = ZERO
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index 91dc705..a4abb66 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -37,6 +37,9 @@
double precision :: tCPU
double precision, external :: wtime
+ ! for Roland_Sylvain integrals
+ double precision :: g_x, g_y, g_z, G_xx, G_yy, G_zz, G_xy, G_xz, G_yz, real_altitude_of_observ_point
+
!--- print number of points and elements in the mesh for each region
if(myrank == 0) then
@@ -49,38 +52,56 @@
! check total Earth mass
if(NCHUNKS == 6) then
write(IMAIN,*)
- write(IMAIN,*) 'computed total Earth mass for this density model and mesh: ',Earth_mass_total
+ write(IMAIN,*) 'computed total Earth mass for this density model and mesh: ',Earth_mass_total,' kg'
write(IMAIN,*) ' (should be not too far from 5.97E+24 kg)'
write(IMAIN,*)
! take into account the fact that dimensions have been non-dimensionalized by dividing them by R_EARTH
- write(IMAIN,*) 'average density for this density model and mesh: ',Earth_mass_total / (volume_total * R_EARTH**3)
+ write(IMAIN,*) 'average density for this density model and mesh: ',Earth_mass_total / (volume_total * R_EARTH**3),' kg/m3'
write(IMAIN,*) ' (should be not too far from 5514 kg/m3)'
endif
!! DK DK for Roland_Sylvain
! Roland_Sylvain integrals
if(ROLAND_SYLVAIN) then
+
+! in m.s-2
+ g_x = Roland_Sylvain_integr_total(1)
+ g_y = Roland_Sylvain_integr_total(2)
+ g_z = Roland_Sylvain_integr_total(3)
+
+! in Eotvos = 1.e+9 s-2
+ G_xx = Roland_Sylvain_integr_total(4) * SI_UNITS_TO_EOTVOS
+ G_yy = Roland_Sylvain_integr_total(5) * SI_UNITS_TO_EOTVOS
+ G_zz = Roland_Sylvain_integr_total(6) * SI_UNITS_TO_EOTVOS
+ G_xy = Roland_Sylvain_integr_total(7) * SI_UNITS_TO_EOTVOS
+ G_xz = Roland_Sylvain_integr_total(8) * SI_UNITS_TO_EOTVOS
+ G_yz = Roland_Sylvain_integr_total(9) * SI_UNITS_TO_EOTVOS
+
write(IMAIN,*)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 1 g_x = ',Roland_Sylvain_integr_total(1)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 2 g_y = ',Roland_Sylvain_integr_total(2)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 3 g_z = ',Roland_Sylvain_integr_total(3)
+ write(IMAIN,*) 'computed total Roland_Sylvain integral g_x = ',g_x,' m.s-2'
+ write(IMAIN,*) 'computed total Roland_Sylvain integral g_y = ',g_y,' m.s-2'
+ write(IMAIN,*) 'computed total Roland_Sylvain integral g_z = ',g_z,' m.s-2'
+ write(IMAIN,*)
+ write(IMAIN,*) 'computed norm of g vector = ',sqrt(g_x**2 + g_y**2 + g_z**2),' m.s-2'
+
+ real_altitude_of_observ_point = sqrt(xobs**2 + yobs**2 + zobs**2)
+! gravity force vector norm decays approximately as (r / r_prime)^2 above the surface of the Earth
+ write(IMAIN,*) ' (should be not too far from ', &
+ sngl(STANDARD_GRAVITY_EARTH * (R_UNIT_SPHERE / real_altitude_of_observ_point)**2),' m.s-2)'
+
write(IMAIN,*)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 4 G_xx = ',Roland_Sylvain_integr_total(4)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 5 G_yy = ',Roland_Sylvain_integr_total(5)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 6 G_zz = ',Roland_Sylvain_integr_total(6)
+ write(IMAIN,*) 'computed total Roland_Sylvain integral G_xx = ',G_xx,' Eotvos'
+ write(IMAIN,*) 'computed total Roland_Sylvain integral G_yy = ',G_yy,' Eotvos'
+ write(IMAIN,*) 'computed total Roland_Sylvain integral G_zz = ',G_zz,' Eotvos'
write(IMAIN,*)
- write(IMAIN,*) 'G tensor should be traceless, G_xx + G_yy + G_zz = 0'
- write(IMAIN,*) 'actual sum obtained = ', &
- Roland_Sylvain_integr_total(4) + Roland_Sylvain_integr_total(5) + Roland_Sylvain_integr_total(6)
- if(max(abs(Roland_Sylvain_integr_total(4)),abs(Roland_Sylvain_integr_total(5)),&
- abs(Roland_Sylvain_integr_total(6))) > TINYVAL) &
- write(IMAIN,*) ' i.e., ',sngl(100.d0*(Roland_Sylvain_integr_total(4) + Roland_Sylvain_integr_total(5) + &
- Roland_Sylvain_integr_total(6)) / max(abs(Roland_Sylvain_integr_total(4)),abs(Roland_Sylvain_integr_total(5)),&
- abs(Roland_Sylvain_integr_total(6)))),'% of max(abs(Gxx),abs(Gyy),abs(Gzz))'
+ write(IMAIN,*) 'G tensor should be traceless, G_xx + G_yy + G_zz = 0.'
+ write(IMAIN,*) 'Actual sum obtained = ',G_xx + G_yy + G_zz
+ if(max(abs(G_xx),abs(G_yy),abs(G_zz)) > TINYVAL) write(IMAIN,*) ' i.e., ', &
+ sngl(100.d0*abs(G_xx + G_yy + G_zz) / max(abs(G_xx),abs(G_yy),abs(G_zz))),'% of max(abs(Gxx),abs(Gyy),abs(Gzz))'
write(IMAIN,*)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 7 G_xy = ',Roland_Sylvain_integr_total(7)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 8 G_xz = ',Roland_Sylvain_integr_total(8)
- write(IMAIN,*) 'computed total Roland_Sylvain integral 9 G_yz = ',Roland_Sylvain_integr_total(9)
+ write(IMAIN,*) 'computed total Roland_Sylvain integral G_xy = ',G_xy,' Eotvos'
+ write(IMAIN,*) 'computed total Roland_Sylvain integral G_xz = ',G_xz,' Eotvos'
+ write(IMAIN,*) 'computed total Roland_Sylvain integral G_yz = ',G_yz,' Eotvos'
endif
More information about the CIG-COMMITS
mailing list