[cig-commits] [commit] devel, master: fixed a significant bug in the calculation of the total volume of the mesh: some central cube elements were counted several times because IFLAG_IN_FICTITIOUS_CUBE elements were not excluded. This volume is computed in the mesher for information purposes only and thus seismograms were not affected. (638b0ef)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:12:31 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 638b0efa6afb9af31fdc4c4c558afbe333304dc2
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed Apr 30 23:44:59 2014 +0200
fixed a significant bug in the calculation of the total volume of the mesh: some central cube elements were counted several times because IFLAG_IN_FICTITIOUS_CUBE elements were not excluded. This volume is computed in the mesher for information purposes only and thus seismograms were not affected.
>---------------------------------------------------------------
638b0efa6afb9af31fdc4c4c558afbe333304dc2
src/meshfem3D/compute_area.f90 | 30 ++----------------------------
src/meshfem3D/compute_volumes.f90 | 31 ++++++++++++++++++++++++++-----
src/meshfem3D/create_regions_mesh.F90 | 27 ++++++++++++++++++---------
src/meshfem3D/finalize_mesher.f90 | 13 ++-----------
4 files changed, 48 insertions(+), 53 deletions(-)
diff --git a/src/meshfem3D/compute_area.f90 b/src/meshfem3D/compute_area.f90
index 4605b4b..61fd234 100644
--- a/src/meshfem3D/compute_area.f90
+++ b/src/meshfem3D/compute_area.f90
@@ -84,7 +84,8 @@
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,*) ' similar area (central cube): ',dble(NCHUNKS)*(2.*(R_CENTRAL_CUBE / R_EARTH)/sqrt(3.))**2
+ 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
@@ -95,30 +96,3 @@
end subroutine compute_area
-!=====================================================================
-
- ! computes total Earth mass
-
- subroutine compute_total_Earth_mass(myrank,Earth_mass_local,Earth_mass_total)
-
- use meshfem3D_models_par
-
- implicit none
-
- integer :: myrank
-
- double precision :: Earth_mass_local
- double precision :: Earth_mass_total
-
- ! local parameters
- double precision :: Earth_mass_total_region
-
- ! use an MPI reduction to compute the total Earth mass
- Earth_mass_total_region = ZERO
- call sum_all_dp(Earth_mass_local,Earth_mass_total_region)
-
- ! sum volume over all the regions
- if(myrank == 0) Earth_mass_total = Earth_mass_total + Earth_mass_total_region
-
- end subroutine compute_total_Earth_mass
-
diff --git a/src/meshfem3D/compute_volumes.f90 b/src/meshfem3D/compute_volumes.f90
index 7610df7..5ebbd0e 100644
--- a/src/meshfem3D/compute_volumes.f90
+++ b/src/meshfem3D/compute_volumes.f90
@@ -28,7 +28,7 @@
subroutine 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)
+ NSPEC2D_BOTTOM,jacobian2D_bottom,NSPEC2D_TOP,jacobian2D_top,idoubling)
use constants
@@ -39,6 +39,8 @@
integer :: nspec
double precision :: wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ)
+ integer,dimension(nspec) :: idoubling
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore
@@ -58,6 +60,10 @@
! calculates volume of all elements in mesh
do ispec = 1,nspec
+
+ ! suppress fictitious elements in central cube
+ if(idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
@@ -110,21 +116,24 @@
!=====================================================================
- ! compute Earth mass of that part of the slice
+ ! compute Earth mass of that part of the slice and then total Earth mass
- subroutine compute_Earth_mass(Earth_mass_local, &
+ subroutine compute_Earth_mass(myrank,Earth_mass_local,Earth_mass_total, &
nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore)
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling)
use constants
implicit none
- double precision :: Earth_mass_local
+ double precision :: Earth_mass_local,Earth_mass_total
+ integer :: myrank
integer :: nspec
double precision :: wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ)
+ integer,dimension(nspec) :: idoubling
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore
@@ -134,6 +143,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
! 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
@@ -143,6 +153,10 @@
! calculates volume of all elements in mesh
do ispec = 1,nspec
+
+ ! suppress fictitious elements in central cube
+ if(idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
+
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
@@ -174,5 +188,12 @@
! 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
+ ! use an MPI reduction to compute the total Earth mass
+ Earth_mass_total_region = ZERO
+ call sum_all_dp(Earth_mass_local,Earth_mass_total_region)
+
+ ! sum volume over all the regions
+ if(myrank == 0) Earth_mass_total = Earth_mass_total + Earth_mass_total_region
+
end subroutine compute_Earth_mass
diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
index 4fd3ba6..1d2a894 100644
--- a/src/meshfem3D/create_regions_mesh.F90
+++ b/src/meshfem3D/create_regions_mesh.F90
@@ -392,7 +392,8 @@
call crm_free_MPI_arrays(iregion_code)
! boundary mesh for MOHO, 400 and 670 discontinuities
- if (SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
+!! DK DK for Roland_Sylvain
+ if (SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE .and. .not. ROLAND_SYLVAIN) then
! user output
call synchronize_all()
if( myrank == 0) then
@@ -413,27 +414,35 @@
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)
+ 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 Earth mass of that part of the slice
- call compute_Earth_mass(Earth_mass_local, &
+ ! 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, &
nspec,wxgll,wygll,wzgll,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore)
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling)
- ! computes total Earth mass
- call compute_total_Earth_mass(myrank,Earth_mass_local,Earth_mass_total)
+!! DK DK for Roland_Sylvain
+ if(ROLAND_SYLVAIN) then
+ call synchronize_all()
+ if( myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) ' ...computing Roland_Sylvain integrals'
+ call flush_IMAIN()
+ endif
+ endif
! create AVS or DX mesh data for the slices
- if(SAVE_MESH_FILES) then
+!! DK DK for Roland_Sylvain
+ if(SAVE_MESH_FILES .and. .not. ROLAND_SYLVAIN) then
! user output
call synchronize_all()
if( myrank == 0) then
write(IMAIN,*)
- write(IMAIN,*) ' ...saving AVS mesh files'
+ write(IMAIN,*) ' ...saving AVS or DX mesh files'
call flush_IMAIN()
endif
call crm_save_mesh_files(nspec,npointot,iregion_code)
diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90
index d717443..9eda5b6 100644
--- a/src/meshfem3D/finalize_mesher.f90
+++ b/src/meshfem3D/finalize_mesher.f90
@@ -43,17 +43,8 @@
! check volume of chunk
write(IMAIN,*)
write(IMAIN,*) 'calculated volume: ',volume_total
- if(.not. TOPOGRAPHY) then
- ! take the central cube into account
- ! it is counted 6 times because of the fictitious elements
- if(INCLUDE_CENTRAL_CUBE) then
- write(IMAIN,*) ' exact volume: ', &
- dble(NCHUNKS)*((4.0d0/3.0d0)*PI*(R_UNIT_SPHERE**3)+5.*(2.*(R_CENTRAL_CUBE/R_EARTH)/sqrt(3.))**3)/6.d0
- else
- write(IMAIN,*) ' exact volume: ', &
- dble(NCHUNKS)*((4.0d0/3.0d0)*PI*(R_UNIT_SPHERE**3)-(2.*(R_CENTRAL_CUBE/R_EARTH)/sqrt(3.))**3)/6.d0
- endif
- endif
+ if(NCHUNKS == 6 .and. .not. ELLIPTICITY .and. .not. TOPOGRAPHY) &
+ write(IMAIN,*) ' exact volume: ',(4.0d0/3.0d0)*PI*(R_UNIT_SPHERE**3)
! check total Earth mass
if(NCHUNKS == 6) then
More information about the CIG-COMMITS
mailing list