Skip to content

Commit ae4c8f9

Browse files
danielpeterdanielpeter
committed
updates PML element assignments (to have at least a single element layer at boundaries for meshes with doubling layers)
1 parent 73af11e commit ae4c8f9

File tree

1 file changed

+239
-56
lines changed

1 file changed

+239
-56
lines changed

‎src/meshfem3D/create_CPML_regions.f90‎

Lines changed: 239 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
3535
! create the different regions of the mesh
3636
use constants, only: IMAIN,CUSTOM_REAL,SMALL_PERCENTAGE_TOLERANCE, &
3737
CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ, &
38-
PI,MAX_STRING_LEN,NDIM,myrank
38+
PI,TINYVAL,HUGEVAL,MAX_STRING_LEN,NDIM,myrank
3939

4040
use constants_meshfem, only: NGLLX_M,NGLLY_M,NGLLZ_M
4141

@@ -50,7 +50,10 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
5050
! local parameters
5151
real(kind=CUSTOM_REAL) :: xmin,xmax,ymin,ymax,zmin,zmax,limit
5252
real(kind=CUSTOM_REAL) :: xmin_all,xmax_all,ymin_all,ymax_all,zmin_all,zmax_all
53-
53+
double precision :: elem_size_x_min,elem_size_x_max
54+
double precision :: elem_size_y_min,elem_size_y_max
55+
double precision :: elem_size_z_min,elem_size_z_max
56+
double precision :: elem_size,elem_size_all
5457
logical, dimension(:), allocatable :: is_X_CPML,is_Y_CPML,is_Z_CPML
5558

5659
integer :: nspec_CPML_total,nspec_total
@@ -59,6 +62,14 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
5962
integer :: ier
6063
character(len=MAX_STRING_LEN) :: filename
6164

65+
! conversion factor from degrees to m (1 degree = 6371.d0 * PI/180 = 111.1949 km) and vice versa
66+
double precision, parameter :: DEGREES_TO_METERS = 6371000.d0 * PI/180.d0
67+
double precision, parameter :: METERS_TO_DEGREES = 1.d0 / (6371000.d0 * PI/180.d0)
68+
69+
! for float comparisons
70+
double precision, parameter :: TOL_EPS = 1.19d-7 ! machine precision for single precision (32-bit) floats
71+
double precision :: rtol ! relative tolerance
72+
6273
! CPML allocation
6374
allocate(is_CPML(nspec),stat=ier)
6475
if (ier /= 0) call exit_MPI_without_rank('error allocating array 1310')
@@ -111,14 +122,14 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
111122
write(IMAIN,*)
112123
endif
113124
! checks
114-
if (THICKNESS_OF_X_PML >= 360.0) stop 'Error invalid degree value for THICKNESS_OF_X_PML'
115-
if (THICKNESS_OF_Y_PML >= 360.0) stop 'Error invalid degree value for THICKNESS_OF_Y_PML'
116-
if (THICKNESS_OF_Z_PML >= 360.0) stop 'Error invalid degree value for THICKNESS_OF_Z_PML'
117-
118-
! converts thickness to m (1 degree = 6371.d0 * PI/180 = 111.1949 km
119-
THICKNESS_OF_X_PML = THICKNESS_OF_X_PML * 6371000.d0 * PI/180.d0
120-
THICKNESS_OF_Y_PML = THICKNESS_OF_Y_PML * 6371000.d0 * PI/180.d0
121-
THICKNESS_OF_Z_PML = THICKNESS_OF_Z_PML * 6371000.d0 * PI/180.d0
125+
if (THICKNESS_OF_X_PML >= 360.d0) stop 'Error invalid degree value for THICKNESS_OF_X_PML'
126+
if (THICKNESS_OF_Y_PML >= 360.d0) stop 'Error invalid degree value for THICKNESS_OF_Y_PML'
127+
if (THICKNESS_OF_Z_PML >= 360.d0) stop 'Error invalid degree value for THICKNESS_OF_Z_PML'
128+
129+
! converts thickness to m
130+
THICKNESS_OF_X_PML = THICKNESS_OF_X_PML * DEGREES_TO_METERS
131+
THICKNESS_OF_Y_PML = THICKNESS_OF_Y_PML * DEGREES_TO_METERS
132+
THICKNESS_OF_Z_PML = THICKNESS_OF_Z_PML * DEGREES_TO_METERS
122133
if (myrank == 0) then
123134
write(IMAIN,*) ' using UTM projection, thickness converted to meters:'
124135
write(IMAIN,*) ' THICKNESS_OF_X_PML (in m) = ',sngl(THICKNESS_OF_X_PML)
@@ -157,6 +168,14 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
157168
is_Y_CPML(:) = .false.
158169
is_Z_CPML(:) = .false.
159170

171+
! stats
172+
elem_size_x_min = HUGEVAL
173+
elem_size_x_max = 0.d0
174+
elem_size_y_min = HUGEVAL
175+
elem_size_y_max = 0.d0
176+
elem_size_z_min = HUGEVAL
177+
elem_size_z_max = 0.d0
178+
160179
do ispec = 1,nspec
161180
! corner points
162181
i1 = ibool(1,1,1,ispec)
@@ -168,64 +187,189 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
168187
i7 = ibool(NGLLX_M,NGLLY_M,NGLLZ_M,ispec)
169188
i8 = ibool(1,NGLLY_M,NGLLZ_M,ispec)
170189

171-
! Xmin CPML
172-
limit = xmin_all + THICKNESS_OF_X_PML * SMALL_PERCENTAGE_TOLERANCE
173-
if ( nodes_coords(i1,1) < limit .and. nodes_coords(i2,1) < limit .and. &
174-
nodes_coords(i3,1) < limit .and. nodes_coords(i4,1) < limit .and. &
175-
nodes_coords(i5,1) < limit .and. nodes_coords(i6,1) < limit .and. &
176-
nodes_coords(i7,1) < limit .and. nodes_coords(i8,1) < limit) then
177-
is_X_CPML(ispec) = .true.
178-
endif
179-
180-
! Xmax CPML
181-
limit = xmax_all - THICKNESS_OF_X_PML * SMALL_PERCENTAGE_TOLERANCE
182-
if ( nodes_coords(i1,1) > limit .and. nodes_coords(i2,1) > limit .and. &
183-
nodes_coords(i3,1) > limit .and. nodes_coords(i4,1) > limit .and. &
184-
nodes_coords(i5,1) > limit .and. nodes_coords(i6,1) > limit .and. &
185-
nodes_coords(i7,1) > limit .and. nodes_coords(i8,1) > limit) then
186-
is_X_CPML(ispec) = .true.
187-
endif
190+
! note: if user THICKNESS_OF_X_PML/THICKNESS_OF_Y_PML/THICKNESS_OF_Z_PML was specified as > 0,
191+
! then we add at least the outer most elements of the mesh, even if these elements are larger
192+
! than the specified thickness.
193+
! this avoids having no PML elements at the mesh boundary when doubling layers are used
194+
! and element sizes vary with depth.
195+
196+
! mesh sides along X
197+
if (THICKNESS_OF_X_PML > 0.d0) then
198+
! Xmin CPML
199+
limit = xmin_all + THICKNESS_OF_X_PML * SMALL_PERCENTAGE_TOLERANCE
200+
if ( nodes_coords(i1,1) < limit .and. nodes_coords(i2,1) < limit .and. &
201+
nodes_coords(i3,1) < limit .and. nodes_coords(i4,1) < limit .and. &
202+
nodes_coords(i5,1) < limit .and. nodes_coords(i6,1) < limit .and. &
203+
nodes_coords(i7,1) < limit .and. nodes_coords(i8,1) < limit) then
204+
is_X_CPML(ispec) = .true.
205+
endif
206+
! Xmax CPML
207+
limit = xmax_all - THICKNESS_OF_X_PML * SMALL_PERCENTAGE_TOLERANCE
208+
if ( nodes_coords(i1,1) > limit .and. nodes_coords(i2,1) > limit .and. &
209+
nodes_coords(i3,1) > limit .and. nodes_coords(i4,1) > limit .and. &
210+
nodes_coords(i5,1) > limit .and. nodes_coords(i6,1) > limit .and. &
211+
nodes_coords(i7,1) > limit .and. nodes_coords(i8,1) > limit) then
212+
is_X_CPML(ispec) = .true.
213+
endif
188214

189-
! Ymin CPML
190-
limit = ymin_all + THICKNESS_OF_Y_PML * SMALL_PERCENTAGE_TOLERANCE
191-
if ( nodes_coords(i1,2) < limit .and. nodes_coords(i2,2) < limit .and. &
192-
nodes_coords(i3,2) < limit .and. nodes_coords(i4,2) < limit .and. &
193-
nodes_coords(i5,2) < limit .and. nodes_coords(i6,2) < limit .and. &
194-
nodes_coords(i7,2) < limit .and. nodes_coords(i8,2) < limit) then
195-
is_Y_CPML(ispec) = .true.
215+
! check if on X-min side (i==1)
216+
if (abs(xmin_all) > 0.0_CUSTOM_REAL) then
217+
rtol = abs(xmin_all) * TOL_EPS ! float comparison tolerance
218+
else
219+
rtol = TOL_EPS ! machine precision
220+
endif
221+
if ( abs(nodes_coords(i1,1) - xmin_all) < rtol &
222+
.or. abs(nodes_coords(i4,1) - xmin_all) < rtol &
223+
.or. abs(nodes_coords(i5,1) - xmin_all) < rtol &
224+
.or. abs(nodes_coords(i8,1) - xmin_all) < rtol ) then
225+
! element has corner(s) on X-min side
226+
is_X_CPML(ispec) = .true.
227+
endif
228+
! check if on X-max side (i==NGLLX)
229+
if (abs(xmax_all) > 0.0_CUSTOM_REAL) then
230+
rtol = abs(xmax_all) * TOL_EPS ! float comparison tolerance
231+
else
232+
rtol = TOL_EPS ! machine precision
233+
endif
234+
if ( abs(nodes_coords(i2,1) - xmax_all) < rtol &
235+
.or. abs(nodes_coords(i3,1) - xmax_all) < rtol &
236+
.or. abs(nodes_coords(i6,1) - xmax_all) < rtol &
237+
.or. abs(nodes_coords(i7,1) - xmax_all) < rtol ) then
238+
! element has corner(s) on X-max side
239+
is_X_CPML(ispec) = .true.
240+
endif
196241
endif
197242

198-
! Ymax CPML
199-
limit = ymax_all - THICKNESS_OF_Y_PML * SMALL_PERCENTAGE_TOLERANCE
200-
if ( nodes_coords(i1,2) > limit .and. nodes_coords(i2,2) > limit .and. &
201-
nodes_coords(i3,2) > limit .and. nodes_coords(i4,2) > limit .and. &
202-
nodes_coords(i5,2) > limit .and. nodes_coords(i6,2) > limit .and. &
203-
nodes_coords(i7,2) > limit .and. nodes_coords(i8,2) > limit) then
204-
is_Y_CPML(ispec) = .true.
205-
endif
243+
! mesh sides along Y
244+
if (THICKNESS_OF_Y_PML > 0.d0) then
245+
! Ymin CPML
246+
limit = ymin_all + THICKNESS_OF_Y_PML * SMALL_PERCENTAGE_TOLERANCE
247+
if ( nodes_coords(i1,2) < limit .and. nodes_coords(i2,2) < limit .and. &
248+
nodes_coords(i3,2) < limit .and. nodes_coords(i4,2) < limit .and. &
249+
nodes_coords(i5,2) < limit .and. nodes_coords(i6,2) < limit .and. &
250+
nodes_coords(i7,2) < limit .and. nodes_coords(i8,2) < limit) then
251+
is_Y_CPML(ispec) = .true.
252+
endif
253+
! Ymax CPML
254+
limit = ymax_all - THICKNESS_OF_Y_PML * SMALL_PERCENTAGE_TOLERANCE
255+
if ( nodes_coords(i1,2) > limit .and. nodes_coords(i2,2) > limit .and. &
256+
nodes_coords(i3,2) > limit .and. nodes_coords(i4,2) > limit .and. &
257+
nodes_coords(i5,2) > limit .and. nodes_coords(i6,2) > limit .and. &
258+
nodes_coords(i7,2) > limit .and. nodes_coords(i8,2) > limit) then
259+
is_Y_CPML(ispec) = .true.
260+
endif
206261

207-
! Zmin CPML
208-
limit = zmin_all + THICKNESS_OF_Z_PML * SMALL_PERCENTAGE_TOLERANCE
209-
if ( nodes_coords(i1,3) < limit .and. nodes_coords(i2,3) < limit .and. &
210-
nodes_coords(i3,3) < limit .and. nodes_coords(i4,3) < limit .and. &
211-
nodes_coords(i5,3) < limit .and. nodes_coords(i6,3) < limit .and. &
212-
nodes_coords(i7,3) < limit .and. nodes_coords(i8,3) < limit) then
213-
is_Z_CPML(ispec) = .true.
262+
! check if on Y-min side (j==1)
263+
if (abs(ymin_all) > 0.0_CUSTOM_REAL) then
264+
rtol = abs(ymin_all) * TOL_EPS ! float comparison tolerance
265+
else
266+
rtol = TOL_EPS ! machine precision
267+
endif
268+
if ( abs(nodes_coords(i1,2) - ymin_all) < rtol &
269+
.or. abs(nodes_coords(i2,2) - ymin_all) < rtol &
270+
.or. abs(nodes_coords(i5,2) - ymin_all) < rtol &
271+
.or. abs(nodes_coords(i6,2) - ymin_all) < rtol ) then
272+
! element has corner(s) on Y-min side
273+
is_Y_CPML(ispec) = .true.
274+
endif
275+
! check if on Y-max side (j==NGLLY)
276+
if (abs(ymax_all) > 0.0_CUSTOM_REAL) then
277+
rtol = abs(ymax_all) * TOL_EPS ! float comparison tolerance
278+
else
279+
rtol = TOL_EPS ! machine precision
280+
endif
281+
if ( abs(nodes_coords(i3,2) - ymax_all) < rtol &
282+
.or. abs(nodes_coords(i4,2) - ymax_all) < rtol &
283+
.or. abs(nodes_coords(i7,2) - ymax_all) < rtol &
284+
.or. abs(nodes_coords(i8,2) - ymax_all) < rtol ) then
285+
! element has corner(s) on Y-max side
286+
is_Y_CPML(ispec) = .true.
287+
endif
214288
endif
215289

216-
if (PML_INSTEAD_OF_FREE_SURFACE) then
217-
! Zmax CPML
218-
limit = zmax_all - THICKNESS_OF_Z_PML * SMALL_PERCENTAGE_TOLERANCE
219-
if ( nodes_coords(i1,3) > limit .and. nodes_coords(i2,3) > limit .and. &
220-
nodes_coords(i3,3) > limit .and. nodes_coords(i4,3) > limit .and. &
221-
nodes_coords(i5,3) > limit .and. nodes_coords(i6,3) > limit .and. &
222-
nodes_coords(i7,3) > limit .and. nodes_coords(i8,3) > limit) then
290+
! mesh sides along Z
291+
if (THICKNESS_OF_Z_PML > 0.d0) then
292+
! Zmin CPML
293+
limit = zmin_all + THICKNESS_OF_Z_PML * SMALL_PERCENTAGE_TOLERANCE
294+
if ( nodes_coords(i1,3) < limit .and. nodes_coords(i2,3) < limit .and. &
295+
nodes_coords(i3,3) < limit .and. nodes_coords(i4,3) < limit .and. &
296+
nodes_coords(i5,3) < limit .and. nodes_coords(i6,3) < limit .and. &
297+
nodes_coords(i7,3) < limit .and. nodes_coords(i8,3) < limit) then
298+
is_Z_CPML(ispec) = .true.
299+
endif
300+
! check if on Z-min side (k==1)
301+
if (abs(zmin_all) > 0.0_CUSTOM_REAL) then
302+
rtol = abs(zmin_all) * TOL_EPS ! float comparison tolerance
303+
else
304+
rtol = TOL_EPS ! machine precision
305+
endif
306+
if ( abs(nodes_coords(i1,3) - zmin_all) < rtol &
307+
.or. abs(nodes_coords(i2,3) - zmin_all) < rtol &
308+
.or. abs(nodes_coords(i3,3) - zmin_all) < rtol &
309+
.or. abs(nodes_coords(i4,3) - zmin_all) < rtol ) then
310+
! element has corner(s) on X-min side
223311
is_Z_CPML(ispec) = .true.
224312
endif
313+
314+
if (PML_INSTEAD_OF_FREE_SURFACE) then
315+
! Zmax CPML
316+
limit = zmax_all - THICKNESS_OF_Z_PML * SMALL_PERCENTAGE_TOLERANCE
317+
if ( nodes_coords(i1,3) > limit .and. nodes_coords(i2,3) > limit .and. &
318+
nodes_coords(i3,3) > limit .and. nodes_coords(i4,3) > limit .and. &
319+
nodes_coords(i5,3) > limit .and. nodes_coords(i6,3) > limit .and. &
320+
nodes_coords(i7,3) > limit .and. nodes_coords(i8,3) > limit) then
321+
is_Z_CPML(ispec) = .true.
322+
endif
323+
! check if on Z-max side (k==NGLLZ)
324+
if (abs(zmax_all) > 0.0_CUSTOM_REAL) then
325+
rtol = abs(zmax_all) * TOL_EPS ! float comparison tolerance
326+
else
327+
rtol = TOL_EPS ! machine precision
328+
endif
329+
if ( abs(nodes_coords(i5,3) - zmax_all) < rtol &
330+
.or. abs(nodes_coords(i6,3) - zmax_all) < rtol &
331+
.or. abs(nodes_coords(i7,3) - zmax_all) < rtol &
332+
.or. abs(nodes_coords(i8,3) - zmax_all) < rtol ) then
333+
! element has corner(s) on Z-max side
334+
is_Z_CPML(ispec) = .true.
335+
endif
336+
endif
225337
endif
226338

339+
! total counter
227340
if (is_X_CPML(ispec) .or. is_Y_CPML(ispec) .or. is_Z_CPML(ispec)) nspec_CPML = nspec_CPML + 1
228341

342+
! stats
343+
if (is_X_CPML(ispec)) then
344+
! element size in X-direction
345+
elem_size = abs(nodes_coords(i1,1)-nodes_coords(i2,1)) ! (1,1,1)-(NGLLX,1,1)
346+
elem_size = max(elem_size,abs(nodes_coords(i4,1)-nodes_coords(i3,1))) ! (1,NGLLY,1)-(NGLLX,NGLLY,1)
347+
elem_size = max(elem_size,abs(nodes_coords(i5,1)-nodes_coords(i6,1))) ! (1,1,NGLLZ)-(NGLLX,1,NGLLZ)
348+
elem_size = max(elem_size,abs(nodes_coords(i8,1)-nodes_coords(i7,1))) ! (1,NGLLY,NGLLZ)-(NGLLX,NGLLY,NGLLZ)
349+
! overall min/max
350+
elem_size_x_min = min(elem_size_x_min,elem_size)
351+
elem_size_x_max = max(elem_size_x_max,elem_size)
352+
endif
353+
if (is_Y_CPML(ispec)) then
354+
! element size in y-direction
355+
elem_size = abs(nodes_coords(i1,2)-nodes_coords(i4,2)) ! (1,1,1)-(1,NGLLY,1)
356+
elem_size = max(elem_size,abs(nodes_coords(i2,2)-nodes_coords(i3,2))) ! (NGLLX,1,1)-(NGLLX,NGLLY,1)
357+
elem_size = max(elem_size,abs(nodes_coords(i5,2)-nodes_coords(i8,2))) ! (1,1,NGLLZ)-(1,NGLLY,NGLLZ)
358+
elem_size = max(elem_size,abs(nodes_coords(i6,2)-nodes_coords(i7,2))) ! (NGLLX,1,NGLLZ)-(NGLLX,NGLLY,NGLLZ)
359+
! overall min/max
360+
elem_size_y_min = min(elem_size_y_min,elem_size)
361+
elem_size_y_max = max(elem_size_y_max,elem_size)
362+
endif
363+
if (is_Z_CPML(ispec)) then
364+
! element size in y-direction
365+
elem_size = abs(nodes_coords(i1,3)-nodes_coords(i5,3)) ! (1,1,1)-(1,1,NGLLZ)
366+
elem_size = max(elem_size,abs(nodes_coords(i2,3)-nodes_coords(i6,3))) ! (NGLLX,1,1)-(NGLLX,1,NGLLZ)
367+
elem_size = max(elem_size,abs(nodes_coords(i4,3)-nodes_coords(i8,3))) ! (1,NGLLY,1)-(1,NGLLY,NGLLZ)
368+
elem_size = max(elem_size,abs(nodes_coords(i3,3)-nodes_coords(i7,3))) ! (NGLLX,NGLLY,1)-(NGLLX,NGLLY,NGLLZ)
369+
! overall min/max
370+
elem_size_z_min = min(elem_size_z_min,elem_size)
371+
elem_size_z_max = max(elem_size_z_max,elem_size)
372+
endif
229373
enddo
230374

231375
! outputs total number of CPML elements
@@ -242,6 +386,45 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords)
242386
write(IMAIN,*) ' Created a total of ',nspec_CPML_total,' unique CPML elements'
243387
write(IMAIN,*) ' (i.e., ',100. * nspec_CPML_total / real(nspec_total),'% of the mesh)'
244388
write(IMAIN,*)
389+
call flush_IMAIN()
390+
endif
391+
392+
! outputs stats of CPML elements
393+
! X side elements
394+
call min_all_all_dp(elem_size_x_min,elem_size_all)
395+
elem_size_x_min = elem_size_all
396+
call max_all_all_dp(elem_size_x_max,elem_size_all)
397+
elem_size_x_max = elem_size_all
398+
! Y side elements
399+
call min_all_all_dp(elem_size_y_min,elem_size_all)
400+
elem_size_y_min = elem_size_all
401+
call max_all_all_dp(elem_size_y_max,elem_size_all)
402+
elem_size_y_max = elem_size_all
403+
! Z side elements
404+
call min_all_all_dp(elem_size_z_min,elem_size_all)
405+
elem_size_z_min = elem_size_all
406+
call max_all_all_dp(elem_size_z_max,elem_size_all)
407+
elem_size_z_max = elem_size_all
408+
409+
if (myrank == 0) then
410+
write(IMAIN,*) ' found PML elements: min/max element size along X sides (in m) = ', &
411+
sngl(elem_size_x_min),'/',sngl(elem_size_x_max)
412+
write(IMAIN,*) ' min/max element size along Y sides (in m) = ', &
413+
sngl(elem_size_y_min),'/',sngl(elem_size_y_max)
414+
write(IMAIN,*) ' min/max element size along Z sides (in m) = ', &
415+
sngl(elem_size_z_min),'/',sngl(elem_size_z_max)
416+
write(IMAIN,*)
417+
! converts to degrees
418+
if (.not. SUPPRESS_UTM_PROJECTION) then
419+
write(IMAIN,'(a,f8.5,a,f8.5)') ' min/max element size along X sides (in degree) = ',&
420+
sngl(elem_size_x_min*METERS_TO_DEGREES),'/',sngl(elem_size_x_max*METERS_TO_DEGREES)
421+
write(IMAIN,'(a,f8.5,a,f8.5)') ' min/max element size along Y sides (in degree) = ',&
422+
sngl(elem_size_y_min*METERS_TO_DEGREES),'/',sngl(elem_size_y_max*METERS_TO_DEGREES)
423+
write(IMAIN,'(a,f8.5,a,f8.5)') ' min/max element size along Z sides (in degree) = ',&
424+
sngl(elem_size_z_min*METERS_TO_DEGREES),'/',sngl(elem_size_z_max*METERS_TO_DEGREES)
425+
write(IMAIN,*)
426+
endif
427+
call flush_IMAIN()
245428
endif
246429

247430
! allocates arrays

0 commit comments

Comments
 (0)