adds check and suggested PML thickness info

Author: danielpeter

src/generate_databases/pml_set_local_dampingcoeff.f90

@@ -316,10 +316,13 @@ subroutine pml_set_local_dampingcoeff(xstore,ystore,zstore)
     write(IMAIN,*)
     ! wavelength
     if (vp_max_all > 0.d0) then
-      write(IMAIN,*) '    resulting dominant wavelength : ',sngl(vp_max_all/f0_FOR_PML),'m'
-      write(IMAIN,*) '               ratio CPML_width_x : ',sngl(CPML_width_x / (vp_max_all/f0_FOR_PML))
-      write(IMAIN,*) '               ratio CPML_width_y : ',sngl(CPML_width_y / (vp_max_all/f0_FOR_PML))
-      write(IMAIN,*) '               ratio CPML_width_z : ',sngl(CPML_width_z / (vp_max_all/f0_FOR_PML))
+      write(IMAIN,*) '    resulting dominant wavelength (L) : ',sngl(vp_max_all/f0_FOR_PML),'m'
+      ! suggested PML width ~ lambda * 0.8  for degree N==4
+      write(IMAIN,*) '          suggested minimum PML width : ',sngl((vp_max_all/f0_FOR_PML) * 0.8d0),'m'
+      write(IMAIN,*) '                 ratio CPML_width_x/L : ',sngl(CPML_width_x / (vp_max_all/f0_FOR_PML)), &
+                     '(suggested minimum ~ 0.5-0.8)'
+      write(IMAIN,*) '                 ratio CPML_width_y/L : ',sngl(CPML_width_y / (vp_max_all/f0_FOR_PML))
+      write(IMAIN,*) '                 ratio CPML_width_z/L : ',sngl(CPML_width_z / (vp_max_all/f0_FOR_PML))
       write(IMAIN,*)
     endif
     call flush_IMAIN()

src/meshfem3D/create_CPML_regions.f90

@@ -609,9 +609,10 @@ subroutine add_CPML_region_as_extra_layers(nglob)
   ! topology of faces of cube (see similar definition in check_mesh_quality.f90)
   integer :: faces_topo(4,6),faces_topo_midpoints(5,6)
   integer :: iface,face_type
+  integer :: elem_mapping(NGNOD)
+  logical :: is_mapped(NGNOD)
   ! MPI Cartesian topology uses W for West (= XI_MIN), E for East (= XI_MAX), S for South (= ETA_MIN), N for North (= ETA_MAX)
   integer, parameter :: W = 1,E = 2,S = 3,N = 4,B = 5,T = 6        ! B==Bottom, T==Top
-  integer :: elem_mapping(NGNOD)
 
   integer :: p1,p2,p3,p4,p5,p6,p7,p8,p9,ia,iglob
   integer :: factor_x,factor_y,factor_z
@@ -795,6 +796,30 @@ subroutine add_CPML_region_as_extra_layers(nglob)
   faces_topo_midpoints(4,T) = 20
   faces_topo_midpoints(5,T) = 26    ! face midpoint
 
+  ! checks face setup
+  do iface = 1,6
+    ! element mapping
+    call get_element_index_mapping(iface,elem_mapping)
+    ! sets node flags to check that mapping is bijective (unique)
+    is_mapped(:) = .false.
+    do ia = 1,NGNOD
+      is_mapped(elem_mapping(ia)) = .true.
+    enddo
+    ! check
+    if (any(is_mapped .eqv. .false.)) then
+      print *,'Error: setup for element mapping is invalid for face type:',iface
+      print *,'       elem_mapping:'
+      do ia = 1,NGNOD
+        print *,'         ',ia,elem_mapping(ia)
+      enddo
+      print *,'       is mapped:'
+      do ia = 1,NGNOD
+        print *,'         ',ia,is_mapped(ia)
+      enddo
+      stop 'Invalid element mapping setup for PML element extension'
+    endif
+  enddo
+
   ! determine element sizes (all thicknesses given in m here)
   SIZE_OF_X_ELEMENT_TO_ADD = 0.d0
   SIZE_OF_Y_ELEMENT_TO_ADD = 0.d0