[cig-commits] [commit] devel: simpler implementation of volume and area calculations in the mesher (merged two routines) (163240c)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Apr 30 15:28:39 PDT 2014
Repository : ssh://geoshell/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/e27d871b067b7bc5c39200e5904d301226b6e5bd...811ae11b80604736d2845c4e5a062755069fc9a6
>---------------------------------------------------------------
commit 163240c5dd1d548938dd75be09361575e362f90f
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu May 1 00:26:50 2014 +0200
simpler implementation of volume and area calculations in the mesher (merged two routines)
>---------------------------------------------------------------
163240c5dd1d548938dd75be09361575e362f90f
src/meshfem3D/compute_area.f90 | 98 ----------------------
...e_volumes.f90 => compute_volumes_and_areas.f90} | 76 +++++++++++++++--
src/meshfem3D/create_meshes.f90 | 3 +-
src/meshfem3D/create_regions_mesh.F90 | 22 ++---
src/meshfem3D/rules.mk | 3 +-
5 files changed, 77 insertions(+), 125 deletions(-)
diff --git a/src/meshfem3D/compute_area.f90 b/src/meshfem3D/compute_area.f90
deleted file mode 100644
index 61fd234..0000000
--- a/src/meshfem3D/compute_area.f90
+++ /dev/null
@@ -1,98 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
-! --------------------------------------------------
-!
-! Main historical authors: Dimitri Komatitsch and Jeroen Tromp
-! Princeton University, USA
-! and CNRS / University of Marseille, France
-! (there are currently many more authors!)
-! (c) Princeton University and CNRS / University of Marseille, April 2014
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine compute_area(myrank,NCHUNKS,iregion_code, &
- area_local_bottom,area_local_top,&
- volume_local,volume_total, &
- RCMB,RICB,R_CENTRAL_CUBE)
-
- use meshfem3D_models_par
-
- implicit none
-
- integer :: myrank,NCHUNKS,iregion_code
-
- double precision :: area_local_bottom,area_local_top,volume_local
- double precision :: volume_total
- double precision :: RCMB,RICB,R_CENTRAL_CUBE
-
- ! local parameters
- double precision :: volume_total_region,area_total_bottom,area_total_top
-
- ! use an MPI reduction to compute the total area and volume
- volume_total_region = ZERO
- area_total_bottom = ZERO
- area_total_top = ZERO
-
- call sum_all_dp(area_local_bottom,area_total_bottom)
- call sum_all_dp(area_local_top,area_total_top)
- call sum_all_dp(volume_local,volume_total_region)
-
- if(myrank == 0) then
- ! sum volume over all the regions
- volume_total = volume_total + volume_total_region
-
- ! check volume of chunk, and bottom and top area
- write(IMAIN,*)
- write(IMAIN,*) ' calculated top area: ',area_total_top
-
- ! compare to exact theoretical value
- if(NCHUNKS == 6 .and. .not. TOPOGRAPHY) then
- select case(iregion_code)
- case(IREGION_CRUST_MANTLE)
- write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*R_UNIT_SPHERE**2
- case(IREGION_OUTER_CORE)
- write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RCMB/R_EARTH)**2
- case(IREGION_INNER_CORE)
- write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RICB/R_EARTH)**2
- case default
- call exit_MPI(myrank,'incorrect region code')
- end select
- endif
-
- write(IMAIN,*) 'calculated bottom area: ',area_total_bottom
-
- ! compare to exact theoretical value
- if(NCHUNKS == 6 .and. .not. TOPOGRAPHY) then
- select case(iregion_code)
- case(IREGION_CRUST_MANTLE)
- write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RCMB/R_EARTH)**2
- case(IREGION_OUTER_CORE)
- write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RICB/R_EARTH)**2
- case(IREGION_INNER_CORE)
- write(IMAIN,*) ' more or less similar area (central cube): ', &
- dble(NCHUNKS)*(2.*(R_CENTRAL_CUBE / R_EARTH)/sqrt(3.))**2
- case default
- call exit_MPI(myrank,'incorrect region code')
- end select
- endif
- call flush_IMAIN()
-
- endif
-
- end subroutine compute_area
-
diff --git a/src/meshfem3D/compute_volumes.f90 b/src/meshfem3D/compute_volumes_and_areas.f90
similarity index 70%
rename from src/meshfem3D/compute_volumes.f90
rename to src/meshfem3D/compute_volumes_and_areas.f90
index 5ebbd0e..adb3132 100644
--- a/src/meshfem3D/compute_volumes.f90
+++ b/src/meshfem3D/compute_volumes_and_areas.f90
@@ -25,20 +25,25 @@
!
!=====================================================================
- subroutine compute_volumes(volume_local,area_local_bottom,area_local_top, &
- nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
+ 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)
+ NSPEC2D_BOTTOM,jacobian2D_bottom,NSPEC2D_TOP,jacobian2D_top,idoubling, &
+ volume_total,RCMB,RICB,R_CENTRAL_CUBE)
use constants
- implicit none
+ use meshfem3D_models_par
- double precision :: volume_local,area_local_bottom,area_local_top
+ implicit none
integer :: nspec
double precision :: wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ)
+ integer :: myrank,NCHUNKS,iregion_code
+
+ double precision :: volume_total
+ double precision :: RCMB,RICB,R_CENTRAL_CUBE
+
integer,dimension(nspec) :: idoubling
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
@@ -49,6 +54,8 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_TOP) :: jacobian2D_top
! local parameters
+ double precision :: volume_local,area_local_bottom,area_local_top
+ double precision :: volume_total_region,area_total_bottom,area_total_top
double precision :: weight
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
integer :: i,j,k,ispec
@@ -112,13 +119,64 @@
enddo
enddo
- end subroutine compute_volumes
+ ! use an MPI reduction to compute the total area and volume
+ volume_total_region = ZERO
+ area_total_bottom = ZERO
+ area_total_top = ZERO
+
+ call sum_all_dp(area_local_bottom,area_total_bottom)
+ call sum_all_dp(area_local_top,area_total_top)
+ call sum_all_dp(volume_local,volume_total_region)
+
+ if(myrank == 0) then
+ ! sum volume over all the regions
+ volume_total = volume_total + volume_total_region
+
+ ! check volume of chunk, and bottom and top area
+ write(IMAIN,*)
+ write(IMAIN,*) ' calculated top area: ',area_total_top
+
+ ! compare to exact theoretical value
+ if(NCHUNKS == 6 .and. .not. TOPOGRAPHY) then
+ select case(iregion_code)
+ case(IREGION_CRUST_MANTLE)
+ write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*R_UNIT_SPHERE**2
+ case(IREGION_OUTER_CORE)
+ write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RCMB/R_EARTH)**2
+ case(IREGION_INNER_CORE)
+ write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RICB/R_EARTH)**2
+ case default
+ call exit_MPI(myrank,'incorrect region code')
+ end select
+ endif
+
+ write(IMAIN,*) 'calculated bottom area: ',area_total_bottom
+
+ ! compare to exact theoretical value
+ if(NCHUNKS == 6 .and. .not. TOPOGRAPHY) then
+ select case(iregion_code)
+ case(IREGION_CRUST_MANTLE)
+ write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RCMB/R_EARTH)**2
+ case(IREGION_OUTER_CORE)
+ write(IMAIN,*) ' exact area: ',dble(NCHUNKS)*(4.0d0/6.0d0)*PI*(RICB/R_EARTH)**2
+ case(IREGION_INNER_CORE)
+ write(IMAIN,*) ' more or less similar area (central cube): ', &
+ dble(NCHUNKS)*(2.*(R_CENTRAL_CUBE / R_EARTH)/sqrt(3.))**2
+ case default
+ call exit_MPI(myrank,'incorrect region code')
+ end select
+ endif
+ call flush_IMAIN()
+
+ endif
+
+ end subroutine compute_volumes_and_areas
!=====================================================================
! compute Earth mass of that part of the slice and then total Earth mass
- subroutine compute_Earth_mass(myrank,Earth_mass_local,Earth_mass_total, &
+ subroutine compute_Earth_mass(myrank,Earth_mass_total, &
nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling)
@@ -126,7 +184,7 @@
implicit none
- double precision :: Earth_mass_local,Earth_mass_total
+ double precision :: Earth_mass_total
integer :: myrank
integer :: nspec
@@ -143,7 +201,7 @@
double precision :: weight
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
integer :: i,j,k,ispec
- double precision :: Earth_mass_total_region
+ double precision :: Earth_mass_local,Earth_mass_total_region
! 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
diff --git a/src/meshfem3D/create_meshes.f90 b/src/meshfem3D/create_meshes.f90
index 89989dc..11665ce 100644
--- a/src/meshfem3D/create_meshes.f90
+++ b/src/meshfem3D/create_meshes.f90
@@ -40,8 +40,9 @@
iproc_xi = iproc_xi_slice(myrank)
iproc_eta = iproc_eta_slice(myrank)
- ! volume of the slice
+ ! volume of the final mesh, and Earth mass computed in the final mesh
volume_total = ZERO
+ Earth_mass_total = ZERO
! make sure everybody is synchronized
call synchronize_all()
diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
index 1d2a894..b7a3e76 100644
--- a/src/meshfem3D/create_regions_mesh.F90
+++ b/src/meshfem3D/create_regions_mesh.F90
@@ -101,10 +101,6 @@
! local parameters
integer :: ier
integer :: nglob
- ! check area and volume of the final mesh
- double precision :: area_local_bottom,area_local_top
- double precision :: volume_local
- double precision :: Earth_mass_local
! user output
if(myrank == 0 ) then
@@ -410,18 +406,14 @@
endif
- ! compute volume, bottom and top area of that part of the slice
- call compute_volumes(volume_local,area_local_bottom,area_local_top, &
- nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
- NSPEC2D_BOTTOM,jacobian2D_bottom,NSPEC2D_TOP,jacobian2D_top,idoubling)
-
- ! computes total area and volume
- call compute_area(myrank,NCHUNKS,iregion_code, area_local_bottom, &
- area_local_top,volume_local,volume_total,RCMB,RICB,R_CENTRAL_CUBE)
+ ! 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, &
+ 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_local,Earth_mass_total, &
+ ! 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, &
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling)
diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk
index d92f9d0..c3aafe5 100644
--- a/src/meshfem3D/rules.mk
+++ b/src/meshfem3D/rules.mk
@@ -44,10 +44,9 @@ meshfem3D_OBJECTS = \
$O/assemble_MPI_scalar_mesh.check.o \
$O/assemble_MPI_vector_mesh.check.o \
$O/calc_jacobian.check.o \
- $O/compute_area.check.o \
$O/compute_coordinates_grid.check.o \
$O/compute_element_properties.check.o \
- $O/compute_volumes.check.o \
+ $O/compute_volumes_and_areas.check.o \
$O/create_addressing.check.o \
$O/create_central_cube.check.o \
$O/create_central_cube_buffers.check.o \
More information about the CIG-COMMITS
mailing list