Merge pull request SPECFEM#1800 from danielpeter/devel

updates submodule path, mesh output and scripts

Author: danielpeter

.gitmodules

@@ -9,4 +9,4 @@
 	url = https://github.com/wjlei1990/pycmt3d
 [submodule "external_libs/ROCm-HIP-CPU"]
 	path = external_libs/ROCm-HIP-CPU
-	url = git@github.com:ROCm/HIP-CPU.git
+	url = https://github.com/ROCm/HIP-CPU.git

EXAMPLES/applications/Mount_StHelens/convert_lonlat2utm.pl

@@ -87,8 +87,14 @@ sub geo2utm {
   #my($PI=pi);
   my($degrad)=pi/180.;
   my($raddeg)=180.0/pi;
-  my($semimaj)=6378206.4;
-  my($semimin)=6356583.8;
+
+  # Clarke 1866
+  #my($semimaj)=6378206.4;
+  #my($semimin)=6356583.8;
+  # WGS84 (World Geodetic System 1984)
+  my($semimaj)=6378137.0;
+  my($semimin)=6356752.314245;
+
   my($scfa)=0.9996;
 
 #! some extracts about UTM:

src/generate_databases/calc_jacobian.f90

@@ -95,7 +95,7 @@ subroutine calc_jacobian(myrank,xix_elem,xiy_elem,xiz_elem, &
           enddo
           write(filename,'(a,i6.6,a)') trim(OUTPUT_FILES)//'/error_proc',myrank,'_element_with_invalid_jacobian'
           call write_VTK_data_points_elem(NGNOD,xelm,yelm,zelm,jacobian,filename)
-          print *,'  written out:',trim(filename)
+          print *,'  written out: ',trim(filename)
           print *,'Please check your mesh...'
           call exit_MPI(myrank,'Error negative or null 3D Jacobian found')
         endif

src/generate_databases/create_regions_mesh.f90

@@ -959,7 +959,7 @@ end subroutine crm_ext_allocate_arrays
 
   subroutine crm_ext_setup_jacobian(nodes_coords_ext_mesh, nnodes_ext_mesh, elmnts_ext_mesh, nelmnts_ext_mesh)
 
-  use constants, only: myrank,NDIM,NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,MAX_STRING_LEN,OUTPUT_FILES
+  use constants, only: myrank,NDIM,NGLLX,NGLLY,NGLLZ,CUSTOM_REAL,MAX_STRING_LEN,OUTPUT_FILES,HUGEVAL,IMAIN
 
   use generate_databases_par, only: NGNOD
 
@@ -987,6 +987,9 @@ subroutine crm_ext_setup_jacobian(nodes_coords_ext_mesh, nnodes_ext_mesh, elmnts
   double precision, dimension(NGNOD) :: xelm,yelm,zelm
   double precision,parameter :: threshold_zero = 1.e-25
 
+  ! stats
+  real(kind=CUSTOM_REAL) :: jacobian_min,jacobian_min_glob,jacobian_max,jacobian_max_glob
+
   ! debug
   logical, parameter :: DEBUG_ELEMENT = .false.
   character(len=MAX_STRING_LEN) :: filename
@@ -1002,6 +1005,9 @@ subroutine crm_ext_setup_jacobian(nodes_coords_ext_mesh, nnodes_ext_mesh, elmnts
   xix_regular = 0.0_CUSTOM_REAL
   jacobian_regular = 0.0_CUSTOM_REAL
 
+  jacobian_min = +HUGEVAL
+  jacobian_max = -HUGEVAL
+
   ! determines regular elements
   do ispec = 1, nspec
     do ia = 1,NGNOD
@@ -1023,6 +1029,10 @@ subroutine crm_ext_setup_jacobian(nodes_coords_ext_mesh, nnodes_ext_mesh, elmnts
                          etaxstore(1,1,1,ispec_irreg),etaystore(1,1,1,ispec_irreg),etazstore(1,1,1,ispec_irreg), &
                          gammaxstore(1,1,1,ispec_irreg),gammaystore(1,1,1,ispec_irreg),gammazstore(1,1,1,ispec_irreg), &
                          jacobianstore(1,1,1,ispec_irreg),xelm,yelm,zelm,dershape3D)
+
+      ! stats
+      if (minval(jacobianstore(:,:,:,ispec_irreg)) < jacobian_min) jacobian_min = minval(jacobianstore(:,:,:,ispec_irreg))
+      if (maxval(jacobianstore(:,:,:,ispec_irreg)) > jacobian_max) jacobian_max = maxval(jacobianstore(:,:,:,ispec_irreg))
     else
       ! sets flag for regular elements
       any_regular_elem = .true.
@@ -1035,7 +1045,7 @@ subroutine crm_ext_setup_jacobian(nodes_coords_ext_mesh, nnodes_ext_mesh, elmnts
       if (myrank == 0 .and. ispec == 1) then
         write(filename,'(a,i6.6,a)') trim(OUTPUT_FILES)//'/proc',myrank,'_debug_element'
         call write_VTK_data_points_elem(NGNOD,xelm,yelm,zelm,dble(jacobianstore(1,1,1,ispec_irreg)),filename)
-        print *,'  written out:',trim(filename)
+        print *,'  written out: ',trim(filename)
       endif
     endif
   enddo
@@ -1100,6 +1110,20 @@ subroutine crm_ext_setup_jacobian(nodes_coords_ext_mesh, nnodes_ext_mesh, elmnts
     ! saves regular values
     xix_regular = xix_reg(1,1,1)
     jacobian_regular  = jacobian_reg(1,1,1)
+
+    ! stats
+    if (jacobian_regular < jacobian_min) jacobian_min = jacobian_regular
+    if (jacobian_regular > jacobian_max) jacobian_max = jacobian_regular
+
+  endif
+
+  ! stats info
+  call min_all_cr(jacobian_min,jacobian_min_glob)
+  call max_all_cr(jacobian_max,jacobian_max_glob)
+  if (myrank == 0) then
+      write(IMAIN,*) '     mesh jacobian: min/max = ',jacobian_min_glob,'/',jacobian_max_glob
+      write(IMAIN,*)
+      call flush_IMAIN()
   endif
 
   end subroutine crm_ext_setup_jacobian

src/generate_databases/save_arrays_solver.F90

@@ -457,8 +457,9 @@ subroutine save_arrays_solver_files()
 
   ! local parameters
   real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: v_tmp
+  real(kind=CUSTOM_REAL), dimension(:), allocatable :: tmp_elem
   integer,dimension(:),allocatable :: v_tmp_i
-  integer :: ier,i,j
+  integer :: ier,i,j,ispec,ispec_irreg
   integer, dimension(:), allocatable :: iglob_tmp
   integer :: inum, num_points
   character(len=MAX_STRING_LEN) :: filename
@@ -687,6 +688,48 @@ subroutine save_arrays_solver_files()
                              xstore_unique,ystore_unique,zstore_unique,ibool, &
                              qkappa_attenuation_store,filename)
 
+
+  ! outputs jacobian
+  if (SAVE_MESH_FILES_ADDITIONAL) then
+    v_tmp(:,:,:,:) = 0.0_CUSTOM_REAL
+    ! minimum jacobian per element
+    allocate(tmp_elem(nspec), stat=ier)
+    if (ier /= 0) stop 'Error allocating temporary tmp_elem array'
+    tmp_elem(:) = 0.0_CUSTOM_REAL
+    do ispec = 1,nspec
+      ispec_irreg = irregular_element_number(ispec)
+      ! irregular_element_number is 0 only if the element is regular
+      if (ispec_irreg /= 0) then
+        ! irregular element
+        v_tmp(:,:,:,ispec) = jacobianstore(:,:,:,ispec_irreg)
+        ! minimum jacobian
+        tmp_elem(ispec) = minval(jacobianstore(:,:,:,ispec_irreg))
+      else
+        ! regular element
+        v_tmp(:,:,:,ispec) = jacobian_regular
+        ! minimum jacobian
+        tmp_elem(ispec) = jacobian_regular
+      endif
+    enddo
+    ! bin file output
+    !open(unit=IOUT,file=prname(1:len_trim(prname))//'jacobian.bin',status='unknown',form='unformatted',iostat=ier)
+    !if (ier /= 0) stop 'error opening file jacobian.bin'
+    !write(IOUT) v_tmp
+    !close(IOUT)
+    ! VTK file output
+    !filename = prname(1:len_trim(prname))//'jacobian'
+    !call write_VTK_data_gll_cr(nspec,nglob_unique, &
+    !                           xstore_unique,ystore_unique,zstore_unique,ibool, &
+    !                           v_tmp,filename)
+    ! VTU file output - minimum jacobian values per element
+    filename = prname(1:len_trim(prname))//'jacobian_min'
+    call write_VTU_data_elem_cr_binary(nspec,nglob_unique, &
+                                       xstore_unique,ystore_unique,zstore_unique,ibool, &
+                                       tmp_elem,filename)
+    ! free memory
+    deallocate(tmp_elem)
+  endif
+
   ! frees temporary array
   deallocate(v_tmp)
 

src/shared/check_mesh_resolution.f90

@@ -98,7 +98,7 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
   ! checks if IMAIN is opened for writing out
   inquire(unit=IMAIN,opened=IMAIN_opened)
 
-! note: this routine can take a long time for large meshes...
+  ! note: this routine can take a long time for large meshes...
   if (myrank == 0 .and. IMAIN_opened) then
     write(IMAIN,*) "Mesh resolution:"
     call flush_IMAIN()
@@ -170,7 +170,7 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
   call max_all_i(NGLOB_AB,NGLOB_AB_global_max)
   call sum_all_i(NGLOB_AB,NGLOB_AB_global_sum)
 
-! outputs infos
+  ! outputs infos
   if (myrank == 0 .and. IMAIN_opened) then
     write(IMAIN,*)
     write(IMAIN,*) '********'
@@ -287,9 +287,6 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
     if (DT_PRESENT) then
       cmax = vel_max * DT / distance_min
       cmax_glob = max(cmax_glob,cmax)
-
-      ! debug: for vtk output
-      if (SAVE_MESH_FILES) tmp1(ispec) = cmax
     endif
 
     ! suggested timestep
@@ -303,7 +300,18 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
     dt_suggested_glob = min(dt_suggested_glob, dt_suggested)
 
     ! debug: for vtk output
-    if (SAVE_MESH_FILES) tmp2(ispec) = pmax
+    if (SAVE_MESH_FILES) then
+      ! Courant/suggested DT number
+      if (DT_PRESENT) then
+        ! Courant number
+        tmp1(ispec) = cmax
+      else
+        ! suggested dt
+        tmp1(ispec) = dt_suggested
+      endif
+      ! minimum period
+      tmp2(ispec) = pmax
+    endif
   enddo
 
   ! Vp velocity
@@ -520,10 +528,15 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
         call flush_IMAIN()
       endif
 
-      ! Courant number
+      ! Courant/suggested DT number
       if (DT_PRESENT) then
+        ! Courant number
         filename = 'res_Courant_number'
         call write_checkmesh_data_hdf5(filename,tmp1)
+      else
+        ! suggested DT
+        filename = 'res_dt_suggested'
+        call write_checkmesh_data_hdf5(filename,tmp1)
       endif
 
       ! minimum period estimate
@@ -541,12 +554,19 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
         call flush_IMAIN()
       endif
 
-      ! Courant number
+      ! Courant/suggested DT number
       if (DT_PRESENT) then
+        ! Courant number
         filename = trim(prname)//'res_Courant_number'
         call write_VTU_data_elem_cr_binary(NSPEC_AB,NGLOB_AB, &
                                            xstore,ystore,zstore,ibool, &
                                            tmp1,filename)
+      else
+        ! suggested DT
+        filename = trim(prname)//'res_dt_suggested'
+        call write_VTU_data_elem_cr_binary(NSPEC_AB,NGLOB_AB, &
+                                           xstore,ystore,zstore,ibool, &
+                                           tmp1,filename)
       endif
       ! minimum period estimate
       filename = trim(prname)//'res_minimum_period'

src/shared/utm_geo.f90

@@ -219,7 +219,7 @@ subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway)
     ! (Universal) Transverse Mercator projection
     !
     ! page 61, eq. 3-21 for M
-    f1 = (1.d0 - e2/4.d0 - 3.d0*e4/64.d0 - 5.d0*e6/256d0)*rlat
+    f1 = (1.d0 - e2/4.d0 - 3.d0*e4/64.d0 - 5.d0*e6/256.d0)*rlat
     f2 = 3.d0*e2/8.d0 + 3.d0*e4/32.d0 + 45.d0*e6/1024.d0
     f2 = f2 * sin(2.d0*rlat)
     ! corrected: using .. + 45 e6 / 1024 instead of .. * 45 e6 / 1024