[cig-commits] [commit] devel, master: added the calculation of the location of the center of mass of the Earth model and mesh (e6f8bd7)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:18:19 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit e6f8bd78d5b43dfc850934421fcfb76b107b2802
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Fri May 23 03:54:11 2014 +0200
added the calculation of the location of the center of mass of the Earth model and mesh
>---------------------------------------------------------------
e6f8bd78d5b43dfc850934421fcfb76b107b2802
src/meshfem3D/compute_volumes_and_areas.F90 | 54 ++++++++++++++++++-----
src/meshfem3D/create_meshes.f90 | 3 ++
src/meshfem3D/create_regions_mesh.F90 | 11 ++---
src/meshfem3D/finalize_mesher.f90 | 8 ++++
src/meshfem3D/meshfem3D_par.f90 | 3 +-
src/specfem3D/compute_forces_crust_mantle_Dev.F90 | 0
src/specfem3D/compute_forces_inner_core_Dev.F90 | 0
src/specfem3D/compute_forces_outer_core_Dev.F90 | 0
src/specfem3D/multiply_arrays_source.f90 | 0
src/specfem3D/update_displacement_Newmark.f90 | 0
10 files changed, 63 insertions(+), 16 deletions(-)
diff --git a/src/meshfem3D/compute_volumes_and_areas.F90 b/src/meshfem3D/compute_volumes_and_areas.F90
index 90460f0..b5ff673 100644
--- a/src/meshfem3D/compute_volumes_and_areas.F90
+++ b/src/meshfem3D/compute_volumes_and_areas.F90
@@ -25,8 +25,8 @@
!
!=====================================================================
- subroutine compute_volumes_and_areas(myrank,NCHUNKS,iregion_code,nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
+ subroutine compute_volumes_and_areas(myrank,NCHUNKS,iregion_code,nspec,wxgll,wygll,wzgll, &
+ xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
NSPEC2D_BOTTOM,jacobian2D_bottom,NSPEC2D_TOP,jacobian2D_top,idoubling, &
volume_total,RCMB,RICB,R_CENTRAL_CUBE)
@@ -177,7 +177,8 @@
! compute Earth mass of that part of the slice and then total Earth mass
subroutine compute_Earth_mass(myrank,Earth_mass_total, &
- nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
+ Earth_center_of_mass_x_total,Earth_center_of_mass_y_total,Earth_center_of_mass_z_total, &
+ nspec,wxgll,wygll,wzgll,xstore,ystore,zstore,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling)
use constants
@@ -185,6 +186,7 @@
implicit none
double precision :: Earth_mass_total
+ double precision :: Earth_center_of_mass_x_total,Earth_center_of_mass_y_total,Earth_center_of_mass_z_total
integer :: myrank
integer :: nspec
@@ -192,6 +194,8 @@
integer,dimension(nspec) :: idoubling
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore
@@ -199,15 +203,24 @@
! local parameters
double precision :: weight
- real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl,rhol
integer :: i,j,k,ispec
- double precision :: Earth_mass_local,Earth_mass_total_region
+
+ double precision :: Earth_mass_local
+ double precision :: Earth_center_of_mass_x_local,Earth_center_of_mass_y_local,Earth_center_of_mass_z_local
+
+ double precision :: Earth_mass_total_region
+ double precision :: Earth_center_of_mass_x_tot_reg,Earth_center_of_mass_y_tot_reg,Earth_center_of_mass_z_tot_reg
! 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
+ double precision, parameter :: non_dimensionalizing_factor1 = RHOAV*R_EARTH**3
+ double precision, parameter :: non_dimensionalizing_factor2 = non_dimensionalizing_factor1 * R_EARTH
! initializes
Earth_mass_local = ZERO
+ Earth_center_of_mass_x_local = ZERO
+ Earth_center_of_mass_y_local = ZERO
+ Earth_center_of_mass_z_local = ZERO
! calculates volume of all elements in mesh
do ispec = 1,nspec
@@ -236,7 +249,13 @@
- xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ xizl*(etaxl*gammayl-etayl*gammaxl))
- Earth_mass_local = Earth_mass_local + dble(jacobianl)*rhostore(i,j,k,ispec)*weight
+ rhol = rhostore(i,j,k,ispec)
+
+ Earth_mass_local = Earth_mass_local + dble(jacobianl)*rhol*weight
+
+ Earth_center_of_mass_x_local = Earth_center_of_mass_x_local + dble(jacobianl)*rhol*xstore(i,j,k,ispec)*weight
+ Earth_center_of_mass_y_local = Earth_center_of_mass_y_local + dble(jacobianl)*rhol*ystore(i,j,k,ispec)*weight
+ Earth_center_of_mass_z_local = Earth_center_of_mass_z_local + dble(jacobianl)*rhol*zstore(i,j,k,ispec)*weight
enddo
enddo
@@ -244,14 +263,29 @@
enddo
! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized
- Earth_mass_local = Earth_mass_local * non_dimensionalizing_factor
+ Earth_mass_local = Earth_mass_local * non_dimensionalizing_factor1
+
+ Earth_center_of_mass_x_local = Earth_center_of_mass_x_local * non_dimensionalizing_factor2
+ Earth_center_of_mass_y_local = Earth_center_of_mass_y_local * non_dimensionalizing_factor2
+ Earth_center_of_mass_z_local = Earth_center_of_mass_z_local * non_dimensionalizing_factor2
! use an MPI reduction to compute the total Earth mass
Earth_mass_total_region = ZERO
+ Earth_center_of_mass_x_tot_reg = ZERO
+ Earth_center_of_mass_y_tot_reg = ZERO
+ Earth_center_of_mass_z_tot_reg = ZERO
call sum_all_dp(Earth_mass_local,Earth_mass_total_region)
+ call sum_all_dp(Earth_center_of_mass_x_local,Earth_center_of_mass_x_tot_reg)
+ call sum_all_dp(Earth_center_of_mass_y_local,Earth_center_of_mass_y_tot_reg)
+ call sum_all_dp(Earth_center_of_mass_z_local,Earth_center_of_mass_z_tot_reg)
- ! sum volume over all the regions
- if(myrank == 0) Earth_mass_total = Earth_mass_total + Earth_mass_total_region
+ ! sum volume over all the regions
+ if(myrank == 0) then
+ Earth_mass_total = Earth_mass_total + Earth_mass_total_region
+ Earth_center_of_mass_x_total = Earth_center_of_mass_x_total + Earth_center_of_mass_x_tot_reg
+ Earth_center_of_mass_y_total = Earth_center_of_mass_y_total + Earth_center_of_mass_y_tot_reg
+ Earth_center_of_mass_z_total = Earth_center_of_mass_z_total + Earth_center_of_mass_z_tot_reg
+ endif
end subroutine compute_Earth_mass
diff --git a/src/meshfem3D/create_meshes.f90 b/src/meshfem3D/create_meshes.f90
index 9c4e4c3..d557d74 100644
--- a/src/meshfem3D/create_meshes.f90
+++ b/src/meshfem3D/create_meshes.f90
@@ -43,6 +43,9 @@
! and Roland_Sylvain integrals
volume_total = ZERO
Earth_mass_total = ZERO
+ Earth_center_of_mass_x_total = ZERO
+ Earth_center_of_mass_y_total = ZERO
+ Earth_center_of_mass_z_total = ZERO
g_x(:,:,:) = ZERO
g_y(:,:,:) = ZERO
diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
index 080a042..89c605d 100644
--- a/src/meshfem3D/create_regions_mesh.F90
+++ b/src/meshfem3D/create_regions_mesh.F90
@@ -47,8 +47,8 @@
use meshfem3D_par,only: &
ibool,idoubling,xstore,ystore,zstore, &
- IMAIN,volume_total,Earth_mass_total,myrank,LOCAL_PATH, &
- IREGION_CRUST_MANTLE,IREGION_OUTER_CORE,IREGION_INNER_CORE, &
+ IMAIN,volume_total,Earth_mass_total,Earth_center_of_mass_x_total,Earth_center_of_mass_y_total,Earth_center_of_mass_z_total, &
+ myrank,LOCAL_PATH,IREGION_CRUST_MANTLE,IREGION_OUTER_CORE,IREGION_INNER_CORE, &
IFLAG_IN_FICTITIOUS_CUBE, &
NCHUNKS,SAVE_MESH_FILES,ABSORBING_CONDITIONS, &
R_CENTRAL_CUBE,RICB,RCMB, &
@@ -419,14 +419,15 @@
endif
! compute volume, bottom and top area of that part of the slice, and then the total
- call compute_volumes_and_areas(myrank,NCHUNKS,iregion_code,nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
+ call compute_volumes_and_areas(myrank,NCHUNKS,iregion_code,nspec,wxgll,wygll,wzgll, &
+ xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
NSPEC2D_BOTTOM,jacobian2D_bottom,NSPEC2D_TOP,jacobian2D_top,idoubling, &
volume_total,RCMB,RICB,R_CENTRAL_CUBE)
! compute Earth mass of that part of the slice, and then total Earth mass
call compute_Earth_mass(myrank,Earth_mass_total, &
- nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
+ Earth_center_of_mass_x_total,Earth_center_of_mass_y_total,Earth_center_of_mass_z_total, &
+ nspec,wxgll,wygll,wzgll,xstore,ystore,zstore,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling)
!! DK DK for Roland_Sylvain
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index 79fde9d..a0fce6c 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -117,6 +117,14 @@
! 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),' kg/m3'
write(IMAIN,*) ' (should be not too far from 5514 kg/m3)'
+ write(IMAIN,*)
+ write(IMAIN,*) 'position of the center of mass of the Earth for this density model and mesh: '
+ write(IMAIN,*) ' x = ',(Earth_center_of_mass_x_total / Earth_mass_total) / 1000.d0,' km'
+ write(IMAIN,*) ' y = ',(Earth_center_of_mass_y_total / Earth_mass_total) / 1000.d0,' km'
+ write(IMAIN,*) ' z = ',(Earth_center_of_mass_z_total / Earth_mass_total) / 1000.d0,' km'
+ write(IMAIN,*) ' distance to center = ',(sqrt(Earth_center_of_mass_x_total**2 + Earth_center_of_mass_y_total**2 + &
+ Earth_center_of_mass_z_total**2) / Earth_mass_total) / 1000.d0,' km'
+ write(IMAIN,*)
endif
!! DK DK for Roland_Sylvain
diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90
index c3dfc76..871af3d 100644
--- a/src/meshfem3D/meshfem3D_par.f90
+++ b/src/meshfem3D/meshfem3D_par.f90
@@ -156,8 +156,9 @@
! check area and volume of the final mesh
double precision :: volume_total
- ! check Earth mass computed in the final mesh
+ ! check Earth mass and Earth center of mass computed in the final mesh
double precision :: Earth_mass_total
+ double precision :: Earth_center_of_mass_x_total,Earth_center_of_mass_y_total,Earth_center_of_mass_z_total
! arrays containing the positions of the observation points in non-dimensionalized value for Roland_Sylvain integrals
! the 1D equivalenced versions are for the FORCE_VECTORIZATION version of the loops
More information about the CIG-COMMITS
mailing list