[cig-commits] [commit] devel, master: added flags SHIFT_TO_THIS_CENTER_OF_MASS and ONLY_COMPUTE_CENTER_OF_MASS for Roland_Sylvain gravity integrals (02ccaa1)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/specfem3d_globe

On branches: devel,master
Link       : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f

>---------------------------------------------------------------

commit 02ccaa1f2b77ee8e1e64df285c54d6a2ce862d56
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date:   Sun May 25 02:48:30 2014 +0200

    added flags SHIFT_TO_THIS_CENTER_OF_MASS and ONLY_COMPUTE_CENTER_OF_MASS for Roland_Sylvain gravity integrals


>---------------------------------------------------------------

02ccaa1f2b77ee8e1e64df285c54d6a2ce862d56
 setup/constants.h.in                        | 13 ++++++++++
 src/meshfem3D/compute_volumes_and_areas.F90 | 38 +++++++++++++++++++++++------
 2 files changed, 43 insertions(+), 8 deletions(-)

diff --git a/setup/constants.h.in b/setup/constants.h.in
index 1d3e3dc..b96cdbf 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -416,6 +416,19 @@
 ! and results for a 3D Earth with topography)
   logical, parameter :: REUSE_EXISTING_OBSERVATION_SURF = .false.
 
+! only compute the center of mass of the 3D model (very cheap), or also compute the gravity integrals (expensive)
+  logical, parameter :: ONLY_COMPUTE_CENTER_OF_MASS = .false.
+
+! we may want to shift the reference frame to a pre-computed center of mass
+  logical, parameter :: SHIFT_TO_THIS_CENTER_OF_MASS = .true.
+
+! position of the center of mass of the Earth for this density model and mesh,
+! but first convert it to meters and then make it non-dimensional
+  double precision, parameter :: x_shift =   0.606220633681674d0 * 1000.d0 / R_EARTH !    km
+  double precision, parameter :: y_shift =   0.433103991863316d0 * 1000.d0 / R_EARTH !    km
+  double precision, parameter :: z_shift =   0.520078637496872d0 * 1000.d0 / R_EARTH !    km
+!    distance to center =   0.908605697566306       km
+
 ! altitude of the observation points in km
   double precision, parameter :: altitude_of_observation_points = 255.d0
 
diff --git a/src/meshfem3D/compute_volumes_and_areas.F90 b/src/meshfem3D/compute_volumes_and_areas.F90
index b5ff673..7dba8ca 100644
--- a/src/meshfem3D/compute_volumes_and_areas.F90
+++ b/src/meshfem3D/compute_volumes_and_areas.F90
@@ -203,6 +203,7 @@
 
   ! local parameters
   double precision :: weight
+  double precision :: x_meshpoint,y_meshpoint,z_meshpoint
   real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl,rhol
   integer :: i,j,k,ispec
 
@@ -253,9 +254,20 @@
 
           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
+          x_meshpoint = xstore(i,j,k,ispec)
+          y_meshpoint = ystore(i,j,k,ispec)
+          z_meshpoint = zstore(i,j,k,ispec)
+
+!! DK DK for Roland_Sylvain gravity calculations we may want to shift the reference frame to a pre-computed center of mass
+          if(SHIFT_TO_THIS_CENTER_OF_MASS) then
+            x_meshpoint = x_meshpoint - x_shift
+            y_meshpoint = y_meshpoint - y_shift
+            z_meshpoint = z_meshpoint - z_shift
+          endif
+
+          Earth_center_of_mass_x_local = Earth_center_of_mass_x_local + dble(jacobianl)*rhol*x_meshpoint*weight
+          Earth_center_of_mass_y_local = Earth_center_of_mass_y_local + dble(jacobianl)*rhol*y_meshpoint*weight
+          Earth_center_of_mass_z_local = Earth_center_of_mass_z_local + dble(jacobianl)*rhol*z_meshpoint*weight
 
         enddo
       enddo
@@ -336,6 +348,9 @@
   integer :: ix,iy,ichunk
 #endif
 
+  ! if we do not want to compute the gravity integrals, only the center of mass (computed before)
+  if(ONLY_COMPUTE_CENTER_OF_MASS) return
+
   ! calculates volume of all elements in mesh
   do ispec = 1,nspec
 
@@ -378,13 +393,20 @@
                         - xiyl*(etaxl*gammazl-etazl*gammaxl) &
                         + xizl*(etaxl*gammayl-etayl*gammaxl))
 
-    if(CHECK_FOR_NEGATIVE_JACOBIANS .and. jacobianl <= ZERO) stop 'error: negative Jacobian found in integral calculation'
+          if(CHECK_FOR_NEGATIVE_JACOBIANS .and. jacobianl <= ZERO) stop 'error: negative Jacobian found in integral calculation'
+
+          x_meshpoint = xstore(i,j,k,ispec)
+          y_meshpoint = ystore(i,j,k,ispec)
+          z_meshpoint = zstore(i,j,k,ispec)
 
-    x_meshpoint = xstore(i,j,k,ispec)
-    y_meshpoint = ystore(i,j,k,ispec)
-    z_meshpoint = zstore(i,j,k,ispec)
+!! DK DK for Roland_Sylvain gravity calculations we may want to shift the reference frame to a pre-computed center of mass
+          if(SHIFT_TO_THIS_CENTER_OF_MASS) then
+            x_meshpoint = x_meshpoint - x_shift
+            y_meshpoint = y_meshpoint - y_shift
+            z_meshpoint = z_meshpoint - z_shift
+          endif
 
-    common_multiplying_factor = jacobianl * weight * rhostore(i,j,k,ispec)
+          common_multiplying_factor = jacobianl * weight * rhostore(i,j,k,ispec)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!!!!!!!!! beginning of loop on all the data to create