renamed local copies of subroutines that are in src/check_mesh_qualit…

…y/check_mesh_quality.f90 to avoid confusion; simplified the calculation of energy in src/specfem3D/compute_energy.f90 to make it faster

src/check_mesh_quality/check_mesh_quality.f90

@@ -165,7 +165,7 @@ program check_mesh_quality
   allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ))
 
 ! compute the derivatives of the 3D shape functions for a 8-node or 27-node element
-  call get_shape3D(dershape3D,xigll,yigll,zigll,NGNOD)
+  call local_version_of_get_shape3D(dershape3D,xigll,yigll,zigll,NGNOD)
 
   allocate(xelm(NGNOD))
   allocate(yelm(NGNOD))
@@ -304,7 +304,7 @@ program check_mesh_quality
 
     if (mod(ispec,100000) == 0) print *,'processed ',ispec,' elements out of ',NSPEC
 
-    call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
+    call local_version_of_create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
                equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax,distmean)
 
 ! store element number in which the edge of minimum or maximum length is located
@@ -398,7 +398,7 @@ program check_mesh_quality
 
 ! loop on all the elements
   do ispec = 1,NSPEC
-    call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
+    call local_version_of_create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
                equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax,distmean)
 
 ! store skewness in histogram
@@ -490,8 +490,9 @@ program check_mesh_quality
 
 ! create script for Gnuplot histogram files
   open(unit=14,file='OUTPUT_FILES/plot_mesh_quality_histograms.gnu',status='unknown')
-  write(14,*) 'set term wxt'
-  write(14,*) '#set term gif'
+  write(14,*) 'set terminal x11'
+  write(14,*) '#set terminal wxt'
+  write(14,*) '#set terminal gif'
 
   if (DISPLAY_HISTOGRAM_DISTMEAN) then
     write(14,*) '#set output "mesh_quality_histogram_meansize.gif"'
@@ -535,7 +536,7 @@ program check_mesh_quality
 ! loop on all the elements
   do ispec = 1,NSPEC
 
-    call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
+    call local_version_of_create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
                equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax,distmean)
 
     if (iformat == 1) then
@@ -581,7 +582,7 @@ program check_mesh_quality
 
   do ispec = ispec_begin,ispec_end
 
-    call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
+    call local_version_of_create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
                equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax,distmean)
 
     if (iformat == 1) then
@@ -615,7 +616,7 @@ program check_mesh_quality
 ! loop on all the elements
   do ispec = ispec_begin,ispec_end
 
-    call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
+    call local_version_of_create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
                equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax,distmean)
 
     if (iformat == 1) then
@@ -654,7 +655,7 @@ end program check_mesh_quality
 
 ! create mesh quality data for a given 3D spectral element
 
-  subroutine create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
+  subroutine local_version_of_create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MAX,delta_t, &
                equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax,distmean)
 
   use constants
@@ -832,15 +833,15 @@ subroutine create_mesh_quality_data_3D(x,y,z,ibool,ispec,NGNOD,NSPEC,NGLOB,VP_MA
    dist4 = sqrt((xelm(4) - xelm(6))**2 + (yelm(4) - yelm(6))**2 + (zelm(4) - zelm(6))**2)
    diagonal_aspect_ratio = max(dist1,dist2,dist3,dist4) / min(dist1,dist2,dist3,dist4)
 
-  end subroutine create_mesh_quality_data_3D
+  end subroutine local_version_of_create_mesh_quality_data_3D
 
 !
 !=====================================================================
 !
 
 ! 3D shape functions for 8-node or 27-node element
 
-  subroutine get_shape3D(dershape3D,xigll,yigll,zigll,NGNOD)
+  subroutine local_version_of_get_shape3D(dershape3D,xigll,yigll,zigll,NGNOD)
 
   use constants
 
@@ -925,7 +926,7 @@ subroutine get_shape3D(dershape3D,xigll,yigll,zigll,NGNOD)
 
           ! note: put further initialization for NGNOD == 27 into subroutine
           !       to avoid compilation errors in case NGNOD == 8
-          call get_shape3D_27(NGNOD,dershape3D,xi,eta,gamma,i,j,k)
+          call local_version_of_get_shape3D_27(NGNOD,dershape3D,xi,eta,gamma,i,j,k)
 
         endif
 
@@ -959,15 +960,15 @@ subroutine get_shape3D(dershape3D,xigll,yigll,zigll,NGNOD)
     enddo
   enddo
 
-  end subroutine get_shape3D
+  end subroutine local_version_of_get_shape3D
 
 !
 !-------------------------------------------------------------------------------------------------
 !
 
 !--- case of a 3D 27-node element
 
-  subroutine get_shape3D_27(NGNOD,dershape3D,xi,eta,gamma,i,j,k)
+  subroutine local_version_of_get_shape3D_27(NGNOD,dershape3D,xi,eta,gamma,i,j,k)
 
   use constants
 
@@ -1106,7 +1107,7 @@ subroutine get_shape3D_27(NGNOD,dershape3D,xi,eta,gamma,i,j,k)
   dershape3D(2,27,i,j,k)=l2xi*l2peta*l2gamma
   dershape3D(3,27,i,j,k)=l2xi*l2eta*l2pgamma
 
-  end subroutine get_shape3D_27
+  end subroutine local_version_of_get_shape3D_27
 
 !
 !=====================================================================

src/decompose_mesh/fault_scotch.f90

@@ -70,7 +70,7 @@ subroutine read_fault_files(localpath_name)
     read(101,*) nbfaults
   else
     nbfaults = 0
-    print *, 'Par_file_faults not found: assume no faults'
+    print *, 'Par_file_faults not found: assuming that there are no faults'
   endif
   close(101)
 

src/generate_databases/fault_generate_databases.f90

@@ -106,7 +106,7 @@ subroutine fault_read_input(prname,myrank)
     read(IIN_PAR,*) nb
     if (myrank == 0) write(IMAIN,*) '  ... reading ', nb,' faults from file DATA/Par_file_faults'
   else
-    if (myrank == 0) write(IMAIN,*) 'File DATA/Par_file_faults not found: assume no faults'
+    if (myrank == 0) write(IMAIN,*) 'File DATA/Par_file_faults not found: assuming that there are no faults'
     close(IIN_PAR)
   endif
 

src/meshfem3D/check_mesh_quality.f90

@@ -308,8 +308,9 @@ subroutine check_mesh_quality(myrank,VP_MAX,NGLOB,NSPEC,x,y,z,ibool, &
     filename = OUTPUT_FILES(1:len_trim(OUTPUT_FILES))//'plot_mesh_quality_histogram.gnu'
 
     open(unit=14,file=trim(filename),status='unknown')
-    write(14,*) 'set term wxt'
-    write(14,*) '#set term gif'
+    write(14,*) 'set terminal x11'
+    write(14,*) '#set terminal wxt'
+    write(14,*) '#set terminal gif'
     write(14,*) '#set output "mesh_quality_histogram.gif"'
     write(14,*)
     write(14,*) 'set xrange [0:1]'

src/specfem3D/compute_energy.f90

@@ -47,13 +47,13 @@ subroutine compute_energy()
   real(kind=CUSTOM_REAL) :: duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
   real(kind=CUSTOM_REAL) :: duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
 
-  real(kind=CUSTOM_REAL) :: epsilon_xx,epsilon_yy,epsilon_zz,epsilon_xy,epsilon_xz,epsilon_yz
-  real(kind=CUSTOM_REAL) :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
-  real(kind=CUSTOM_REAL) :: vx,vy,vz,pressure
+  double precision :: epsilon_xx,epsilon_yy,epsilon_zz,epsilon_xy,epsilon_xz,epsilon_yz
+  double precision :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+  double precision :: vx,vy,vz,pressure
 
   real(kind=CUSTOM_REAL) :: hp1,hp2,hp3
 
-  real(kind=CUSTOM_REAL) :: lambdal,mul,lambdalplus2mul,rhol,cpl
+  real(kind=CUSTOM_REAL) :: lambdal,mul,lambdalplus2mul,rhol
   real(kind=CUSTOM_REAL) :: kappal
 
   real(kind=CUSTOM_REAL) :: integration_weight
@@ -166,9 +166,9 @@ subroutine compute_energy()
             epsilon_xx = duxdxl
             epsilon_yy = duydyl
             epsilon_zz = duzdzl
-            epsilon_xy = 0.5 * duxdyl_plus_duydxl
-            epsilon_xz = 0.5 * duzdxl_plus_duxdzl
-            epsilon_yz = 0.5 * duzdyl_plus_duydzl
+            epsilon_xy = 0.5d0 * duxdyl_plus_duydxl
+            epsilon_xz = 0.5d0 * duzdxl_plus_duxdzl
+            epsilon_yz = 0.5d0 * duzdyl_plus_duydzl
 
             kappal = kappastore(i,j,k,ispec)
             mul = mustore(i,j,k,ispec)
@@ -190,14 +190,13 @@ subroutine compute_energy()
             integration_weight = wxgll(i)*wygll(j)*wzgll(k)*jacobianl
 
             ! compute kinetic energy  1/2 rho ||v||^2
-            kinetic_energy = kinetic_energy + 0.5d0 * integration_weight * rhol*(veloc(1,iglob)**2 + &
+            kinetic_energy = kinetic_energy + integration_weight * rhol*(veloc(1,iglob)**2 + &
                                  veloc(2,iglob)**2 + veloc(3,iglob)**2)
 
             ! compute potential energy 1/2 sigma_ij epsilon_ij
-            potential_energy = potential_energy + integration_weight * 0.5d0 * &
-              (sigma_xx*epsilon_xx + sigma_xy*epsilon_xy + sigma_xz*epsilon_xz + &
-               sigma_xy*epsilon_xy + sigma_yy*epsilon_yy + sigma_yz*epsilon_yz + &
-               sigma_xz*epsilon_xz + sigma_yz*epsilon_yz + sigma_zz*epsilon_zz)
+            potential_energy = potential_energy + integration_weight * &
+              (sigma_xx*epsilon_xx + sigma_yy*epsilon_yy + sigma_zz*epsilon_zz + &
+               2.d0 * (sigma_xy*epsilon_xy + sigma_xz*epsilon_xz + sigma_yz*epsilon_yz))
 
           enddo
         enddo
@@ -268,7 +267,6 @@ subroutine compute_energy()
 
             rhol = rhostore(i,j,k,ispec)
             kappal = kappastore(i,j,k,ispec)
-            cpl = sqrt(kappal / rhol)
 
             ! velocity is v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi)
             vx = duxdxl / rhol
@@ -281,10 +279,10 @@ subroutine compute_energy()
             integration_weight = wxgll(i)*wygll(j)*wzgll(k)*jacobianl
 
             ! compute kinetic energy  1/2 rho ||v||^2
-            kinetic_energy = kinetic_energy + integration_weight * 0.5d0 * rhol*(vx**2 + vy**2 + vz**2)
+            kinetic_energy = kinetic_energy + integration_weight * rhol*(vx**2 + vy**2 + vz**2)
 
             ! compute potential energy 1/2 sigma_ij epsilon_ij
-            potential_energy = potential_energy + integration_weight * pressure**2 / (2. * rhol * cpl**2)
+            potential_energy = potential_energy + integration_weight * pressure**2 / kappal
 
           enddo
         enddo
@@ -299,6 +297,8 @@ subroutine compute_energy()
   enddo
 
 ! compute the total using a reduction between all the processors
+  kinetic_energy = 0.5d0 * kinetic_energy
+  potential_energy = 0.5d0 * potential_energy
   call sum_all_dp(kinetic_energy,kinetic_energy_glob)
   call sum_all_dp(potential_energy,potential_energy_glob)
   total_energy_glob = kinetic_energy_glob + potential_energy_glob

src/specfem3D/compute_forces_acoustic_NGLL5_fast.f90

@@ -45,7 +45,8 @@ subroutine compute_forces_acoustic_NGLL5_fast(iphase,NSPEC_AB,NGLOB_AB, &
   use specfem_par, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,chi_elem,temp1,temp2,temp3,temp4, &
                          PML_dpotential_dxl,PML_dpotential_dyl,PML_dpotential_dzl, &
                          PML_dpotential_dxl_old,PML_dpotential_dyl_old,PML_dpotential_dzl_old, &
-                         PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new
+                         PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new, &
+                         OUTPUT_ENERGY,USE_TRICK_FOR_BETTER_PRESSURE,myrank
 
   use pml_par, only: is_CPML, spec_to_CPML, potential_dot_dot_acoustic_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl, &
                      rmemory_dpotential_dzl,rmemory_potential_acoustic, &
@@ -103,6 +104,9 @@ subroutine compute_forces_acoustic_NGLL5_fast(iphase,NSPEC_AB,NGLOB_AB, &
     num_elements = nspec_inner_acoustic
   endif
 
+  if (OUTPUT_ENERGY .and. USE_TRICK_FOR_BETTER_PRESSURE) &
+    call exit_mpi(myrank,'USE_TRICK_FOR_BETTER_PRESSURE not implemented for OUTPUT_ENERGY yet, please turn one of them off')
+
   ! loop over spectral elements
   do ispec_p = 1,num_elements
 
@@ -260,7 +264,7 @@ subroutine compute_forces_acoustic_NGLL5_fast(iphase,NSPEC_AB,NGLOB_AB, &
           PML_dpotential_dyl_new(i,j,k) = xiyl*temp1l_new + etayl*temp2l_new + gammayl*temp3l_new
           PML_dpotential_dzl_new(i,j,k) = xizl*temp1l_new + etazl*temp2l_new + gammazl*temp3l_new
 
-          ! density (reciproc)
+          ! reciprocal of density
           rho_invl = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec)
 
           ! for acoustic medium

src/specfem3D/compute_forces_acoustic_NGLLnot5_generic_slow.f90

@@ -45,7 +45,8 @@ subroutine compute_forces_acoustic_NGLLnot5_generic_slow(iphase,NSPEC_AB,NGLOB_A
   use specfem_par, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,chi_elem,temp1,temp2,temp3,temp4, &
                          PML_dpotential_dxl,PML_dpotential_dyl,PML_dpotential_dzl, &
                          PML_dpotential_dxl_old,PML_dpotential_dyl_old,PML_dpotential_dzl_old, &
-                         PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new
+                         PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new, &
+                         OUTPUT_ENERGY,USE_TRICK_FOR_BETTER_PRESSURE,myrank
 
   use pml_par, only: is_CPML, spec_to_CPML, potential_dot_dot_acoustic_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl, &
                      rmemory_dpotential_dzl,rmemory_potential_acoustic, &
@@ -103,6 +104,9 @@ subroutine compute_forces_acoustic_NGLLnot5_generic_slow(iphase,NSPEC_AB,NGLOB_A
     num_elements = nspec_inner_acoustic
   endif
 
+  if (OUTPUT_ENERGY .and. USE_TRICK_FOR_BETTER_PRESSURE) &
+    call exit_mpi(myrank,'USE_TRICK_FOR_BETTER_PRESSURE not implemented for OUTPUT_ENERGY yet, please turn one of them off')
+
   ! loop over spectral elements
   do ispec_p = 1,num_elements
 
@@ -211,7 +215,7 @@ subroutine compute_forces_acoustic_NGLLnot5_generic_slow(iphase,NSPEC_AB,NGLOB_A
           PML_dpotential_dyl_new(i,j,k) = xiyl*temp1l_new + etayl*temp2l_new + gammayl*temp3l_new
           PML_dpotential_dzl_new(i,j,k) = xizl*temp1l_new + etazl*temp2l_new + gammazl*temp3l_new
 
-          ! density (reciproc)
+          ! reciprocal of density
           rho_invl = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec)
 
           ! for acoustic medium

src/specfem3D/compute_forces_poroelastic_calling_routine.f90

@@ -38,6 +38,8 @@ subroutine compute_forces_poroelastic_calling()
 
   integer:: iphase
 
+  if (OUTPUT_ENERGY) call exit_MPI(myrank,'calculation of energy currently not implemented for poroelastic media')
+
 ! distinguishes two runs: for elements in contact with MPI interfaces, and elements within the partitions
   do iphase = 1,2
 

src/specfem3D/iterate_time.F90

@@ -53,19 +53,20 @@ subroutine iterate_time()
   !----  create a Gnuplot script to display the energy curve in log scale
   if (OUTPUT_ENERGY .and. myrank == 0) then
     open(unit=IOUT_ENERGY,file=trim(OUTPUT_FILES)//'plot_energy.gnu',status='unknown',action='write')
-    write(IOUT_ENERGY,*) 'set term wxt'
-    write(IOUT_ENERGY,*) '#set term postscript landscape color solid "Helvetica" 22'
+    write(IOUT_ENERGY,*) 'set terminal x11'
+    write(IOUT_ENERGY,*) '#set terminal wxt'
+    write(IOUT_ENERGY,*) '#set terminal postscript landscape color solid "Helvetica" 22'
     write(IOUT_ENERGY,*) '#set output "energy.ps"'
     write(IOUT_ENERGY,*) 'set logscale y'
     write(IOUT_ENERGY,*) 'set xlabel "Time step number"'
     write(IOUT_ENERGY,*) 'set ylabel "Energy (J)"'
-    write(IOUT_ENERGY,'(a152)') '#plot "energy.dat" us 1:2 t "Kinetic Energy" w l lc 1, "energy.dat" us 1:3 &
+    write(IOUT_ENERGY,'(a152)') 'plot "energy.dat" us 1:2 t "Kinetic Energy" w l lc 1, "energy.dat" us 1:3 &
                          &t "Potential Energy" w l lc 2, "energy.dat" us 1:4 t "Total Energy" w l lc 4'
-    write(IOUT_ENERGY,*) '#pause -1 "Hit any key..."'
-    write(IOUT_ENERGY,*) '#plot "energy.dat" us 1:2 t "Kinetic Energy" w l lc 1'
-    write(IOUT_ENERGY,*) '#pause -1 "Hit any key..."'
-    write(IOUT_ENERGY,*) '#plot "energy.dat" us 1:3 t "Potential Energy" w l lc 2'
-    write(IOUT_ENERGY,*) '#pause -1 "Hit any key..."'
+    write(IOUT_ENERGY,*) 'pause -1 "Hit any key..."'
+    write(IOUT_ENERGY,*) 'plot "energy.dat" us 1:2 t "Kinetic Energy" w l lc 1'
+    write(IOUT_ENERGY,*) 'pause -1 "Hit any key..."'
+    write(IOUT_ENERGY,*) 'plot "energy.dat" us 1:3 t "Potential Energy" w l lc 2'
+    write(IOUT_ENERGY,*) 'pause -1 "Hit any key..."'
     write(IOUT_ENERGY,*) 'plot "energy.dat" us 1:4 t "Total Energy" w l lc 4'
     write(IOUT_ENERGY,*) 'pause -1 "Hit any key..."'
     close(IOUT_ENERGY)

utils/location_of_Matthieu_s_tool_to_validate_sensitivity_kernels_feb2017.txt

@@ -0,0 +1,45 @@
+
+Subject: Tests for SPECFEM3D_GLOBE and ASDF
+From: Matthieu Lefebvre
+Date: 02/21/2017 05:04 PM
+
+Dear all,
+
+Below is the link to the comparison tool.
+
+Best,
+Matthieu
+
+----
+On 2/6/17, 4:43 PM, Matthieu Lefebvre wrote:
+
+    Dear all,
+
+    1. It is half accurate. James did ran allinea map on specfem3d_globe and has shown parallel-hdf5 to produce a lot of network traffic. As a result, it was marked as an item for IBM's CoE team to investigate. They were later confused by the remarks about the optimization of the receivers' location.
+    2. I have shared a utility to validate kernels with Malte. It should be integrated in a Jenkin test very soon. https://github.com/mpbl/CompareAdiosKernels
+    3. See (2). I highly doubt with can use Titan as a CI machine. The way to go is to test independent pieces on relevant but small cases (e.g. regional Greece). We might then define integration tests involving multiple parts of the workflow.
+
+    Best regards,
+    Matthieu
+
+----
+    On 2/3/17, 1:08 PM, Youyi Ruan wrote:
+
+        Dear All,
+
+        Ebru recently reported almost triple runtime (9 min to 25 min with same resolution and length) than usual using latest version of specfem3d_globe and ASDF on Titian. In regard to her reports and to avoid such issue in future, I am writing to collect updates from the whole group about the current existing tests and the tests we need in the following list upon Jeroen's request.
+
+        1, ASDF tests for both correctness and performance. For I/O performance I remember Matthieu have reported that the overhead cost of thousand receivers should within 1 min or so on Titan.   Please let us know if the ASDF tests are exist in current buildbot test otherwise we will need to add the reference ASDF file for test, and maybe timing module to count the I/O cost.   James, Matthieu and David, would you please respond to or comment on this?
+
+        2, ADIOS test. (This part is beyond my knowledge and I will let Matthieu or Dimitri to respond what tests have been implemented and what we need)
+
+        3, Gradient/kernel tests for both correctness and performance. We have discussed this before, it seems a challenge job for buildbot on Tiger, but at least we will need to setup a test on Titan for any major updates.  Matthieu has been working on the testing code comparing two gradients and discussed how to set the threshold on Thursday's group meeting.  Matthieu, would you please give us an update? Thanks.
+
+        Please also comment if you think there is anything else you would suggest for test.
+
+        Best,
+
+        Youyi
+
+
+