[cig-commits] [commit] devel, master: faster way of computing the Roland_Sylvain integrals by pre-computing many common terms (e347ecc)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:13:14 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit e347eccad0cfbea066ff5efc78c06551c29dfebe
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu May 1 17:53:18 2014 +0200
faster way of computing the Roland_Sylvain integrals by pre-computing many common terms
>---------------------------------------------------------------
e347eccad0cfbea066ff5efc78c06551c29dfebe
setup/constants.h.in | 19 ++++----
src/meshfem3D/compute_volumes_and_areas.f90 | 75 ++++++++++++++++++++---------
src/meshfem3D/finalize_mesher.f90 | 2 +-
3 files changed, 63 insertions(+), 33 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 79fa1b9..8a9d6a9 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -413,20 +413,19 @@
! 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
+! 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
!! 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
+! 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
!! 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 :: 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
!------------------------------------------------------
!----------- do not modify anything below -------------
diff --git a/src/meshfem3D/compute_volumes_and_areas.f90 b/src/meshfem3D/compute_volumes_and_areas.f90
index 050408c..ef05c94 100644
--- a/src/meshfem3D/compute_volumes_and_areas.f90
+++ b/src/meshfem3D/compute_volumes_and_areas.f90
@@ -288,10 +288,16 @@
double precision :: jacobianl
integer :: i,j,k,ispec
double precision :: xval,yval,zval
+ double precision :: xval_squared,yval_squared,zval_squared
+ double precision :: x_meshpoint,y_meshpoint,z_meshpoint,rho_meshpoint
double precision :: distance,distance_squared,distance_cubed,distance_fifth_power, &
- one_over_distance_squared,one_over_distance_cubed,one_over_distance_fifth_power
+ three_over_distance_squared,one_over_distance_cubed,three_over_distance_fifth_power
+ double precision :: common_multiplying_factor
- double precision, dimension(9) :: Roland_Sylvain_int_local,Roland_Sylvain_int_total_region,elemental_contribution
+ double precision, dimension(9) :: Roland_Sylvain_int_local,Roland_Sylvain_int_total_region
+ double precision :: elemental_contribution_1,elemental_contribution_2,elemental_contribution_3, &
+ elemental_contribution_4,elemental_contribution_5,elemental_contribution_6, &
+ elemental_contribution_7,elemental_contribution_8,elemental_contribution_9
! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized
! for the gravity vector force, a distance is involved in the dimensions
@@ -333,50 +339,75 @@
- xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ xizl*(etaxl*gammayl-etayl*gammaxl))
- xval = xstore(i,j,k,ispec) - xobs
- yval = ystore(i,j,k,ispec) - yobs
- zval = zstore(i,j,k,ispec) - zobs
+ x_meshpoint = xstore(i,j,k,ispec)
+ y_meshpoint = ystore(i,j,k,ispec)
+ z_meshpoint = zstore(i,j,k,ispec)
-!! DK DK see later if the non-dimensional values are correctly handled here (I will check more carefully)
- distance_squared = xval**2 + yval**2 + zval**2
+ rho_meshpoint = rhostore(i,j,k,ispec)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!! beginning of loop on all the data to create
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ xval = x_meshpoint - x_observation
+ yval = y_meshpoint - y_observation
+ zval = z_meshpoint - z_observation
+
+ xval_squared = xval**2
+ yval_squared = yval**2
+ zval_squared = zval**2
+
+ distance_squared = xval_squared + yval_squared + zval_squared
distance = sqrt(distance_squared)
distance_cubed = distance_squared*distance
distance_fifth_power = distance_squared*distance_cubed
- one_over_distance_squared = 1.d0 / distance_squared
+ three_over_distance_squared = 3.d0 / distance_squared
one_over_distance_cubed = 1.d0 / distance_cubed
-!!!!!!!!!!!!! one_over_distance_fifth_power = 1.d0 / distance_fifth_power
-!! DK DK this will be faster
- one_over_distance_fifth_power = one_over_distance_squared * one_over_distance_cubed
+ three_over_distance_fifth_power = three_over_distance_squared * one_over_distance_cubed
! g_x
- elemental_contribution(1) = xval * one_over_distance_cubed
+ elemental_contribution_1 = xval * one_over_distance_cubed
! g_y
- elemental_contribution(2) = yval * one_over_distance_cubed
+ elemental_contribution_2 = yval * one_over_distance_cubed
! g_z
- elemental_contribution(3) = zval * one_over_distance_cubed
+ elemental_contribution_3 = zval * one_over_distance_cubed
! G_xx
- elemental_contribution(4) = (3.d0 * xval**2 * one_over_distance_squared - 1.d0) * one_over_distance_cubed
+ elemental_contribution_4 = (xval_squared * three_over_distance_squared - 1.d0) * one_over_distance_cubed
! G_yy
- elemental_contribution(5) = (3.d0 * yval**2 * one_over_distance_squared - 1.d0) * one_over_distance_cubed
+ elemental_contribution_5 = (yval_squared * three_over_distance_squared - 1.d0) * one_over_distance_cubed
! G_zz
- elemental_contribution(6) = (3.d0 * zval**2 * one_over_distance_squared - 1.d0) * one_over_distance_cubed
+ elemental_contribution_6 = (zval_squared * three_over_distance_squared - 1.d0) * one_over_distance_cubed
! G_xy
- elemental_contribution(7) = 3.d0 * xval*yval * one_over_distance_fifth_power
+ elemental_contribution_7 = xval*yval * three_over_distance_fifth_power
! G_xz
- elemental_contribution(8) = 3.d0 * xval*zval * one_over_distance_fifth_power
+ elemental_contribution_8 = xval*zval * three_over_distance_fifth_power
! G_yz
- elemental_contribution(9) = 3.d0 * yval*zval * one_over_distance_fifth_power
-
- Roland_Sylvain_int_local(:) = Roland_Sylvain_int_local(:) + jacobianl*weight*rhostore(i,j,k,ispec)*elemental_contribution(:)
+ elemental_contribution_9 = yval*zval * three_over_distance_fifth_power
+
+ common_multiplying_factor = jacobianl * weight * rho_meshpoint
+
+ Roland_Sylvain_int_local(1) = Roland_Sylvain_int_local(1) + common_multiplying_factor*elemental_contribution_1
+ Roland_Sylvain_int_local(2) = Roland_Sylvain_int_local(2) + common_multiplying_factor*elemental_contribution_2
+ Roland_Sylvain_int_local(3) = Roland_Sylvain_int_local(3) + common_multiplying_factor*elemental_contribution_3
+ Roland_Sylvain_int_local(4) = Roland_Sylvain_int_local(4) + common_multiplying_factor*elemental_contribution_4
+ Roland_Sylvain_int_local(5) = Roland_Sylvain_int_local(5) + common_multiplying_factor*elemental_contribution_5
+ Roland_Sylvain_int_local(6) = Roland_Sylvain_int_local(6) + common_multiplying_factor*elemental_contribution_6
+ Roland_Sylvain_int_local(7) = Roland_Sylvain_int_local(7) + common_multiplying_factor*elemental_contribution_7
+ Roland_Sylvain_int_local(8) = Roland_Sylvain_int_local(8) + common_multiplying_factor*elemental_contribution_8
+ Roland_Sylvain_int_local(9) = Roland_Sylvain_int_local(9) + common_multiplying_factor*elemental_contribution_9
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!! end of loop on all the data to create
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
enddo
enddo
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index a4abb66..7fc4f47 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -84,7 +84,7 @@
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)
+ real_altitude_of_observ_point = sqrt(x_observation**2 + y_observation**2 + z_observation**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)'
More information about the CIG-COMMITS
mailing list