[cig-commits] r21765 - seismo/3D/SPECFEM3D/trunk/utils/CPML
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Apr 8 10:06:14 PDT 2013
Author: dkomati1
Date: 2013-04-08 10:06:13 -0700 (Mon, 08 Apr 2013)
New Revision: 21765
Added:
seismo/3D/SPECFEM3D/trunk/utils/CPML/Makefile
Modified:
seismo/3D/SPECFEM3D/trunk/utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers.f90
Log:
done cleaning and making utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers.f90 more general
Added: seismo/3D/SPECFEM3D/trunk/utils/CPML/Makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/CPML/Makefile (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/CPML/Makefile 2013-04-08 17:06:13 UTC (rev 21765)
@@ -0,0 +1,25 @@
+
+# Fortran90 compiler
+F90 = gfortran
+#F90 = ifort
+
+FFLAGS = -O3
+
+default: xconvert_external_layers_to_CPML cleanobj
+
+all: clean default
+
+xconvert_external_layers_to_CPML: convert_external_layers_of_a_given_mesh_to_CPML_layers.o
+ ${F90} $(FFLAGS) -o xconvert_external_layers_to_CPML convert_external_layers_of_a_given_mesh_to_CPML_layers.o
+
+clean:
+ rm -rf *.o xconvert_external_layers_to_CPML
+
+cleanobj:
+ rm -rf *.o
+
+# OBJECTS:
+
+convert_external_layers_of_a_given_mesh_to_CPML_layers.o: convert_external_layers_of_a_given_mesh_to_CPML_layers.f90
+ ${F90} $(FFLAGS) -c -o convert_external_layers_of_a_given_mesh_to_CPML_layers.o convert_external_layers_of_a_given_mesh_to_CPML_layers.f90
+
Modified: seismo/3D/SPECFEM3D/trunk/utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers.f90 2013-04-08 16:54:30 UTC (rev 21764)
+++ seismo/3D/SPECFEM3D/trunk/utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers.f90 2013-04-08 17:06:13 UTC (rev 21765)
@@ -14,7 +14,7 @@
integer :: nspec,npoin
integer :: ispec,ipoin
- integer :: ipoin_read,ispec_read
+ integer :: ipoin_read,ispec_loop,iflag
integer :: i1,i2,i3,i4,i5,i6,i7,i8,number_of_CPML_elements
real, dimension(:), allocatable :: x,y,z
@@ -25,17 +25,46 @@
real :: xread,yread,zread,xmin,xmax,ymin,ymax,zmin,zmax,limit
- logical, parameter :: ALSO_ADD_ON_THE_TOP_SURFACE = .false.
+ logical :: ALSO_ADD_ON_THE_TOP_SURFACE
- real, parameter :: THICKNESS_OF_X_PML = 3722.219 * 3 !! DK DK because we want 3 CPML elements
- real, parameter :: THICKNESS_OF_Y_PML = THICKNESS_OF_X_PML
- real, parameter :: THICKNESS_OF_Z_PML = 3750. * 3 !! DK DK because we want 3 CPML elements
+ real :: THICKNESS_OF_X_PML, THICKNESS_OF_Y_PML, THICKNESS_OF_Z_PML
- real, parameter :: SMALL_PERCENTAGE_TOLERANCE = 1.03 !! to make sure coordinate roundoff problems do not occur
+! to make sure coordinate roundoff problems do not occur, use a tolerance of 1.5%
+ real, parameter :: SMALL_PERCENTAGE_TOLERANCE = 1.015
-!! DK DK dire dans manuel et ici en comment: it is your responsibility to have flat mesh elements of the right thickness aligned
-!! with the coordinate axes in your input CUBIT (or similar) mesh; this code does NOT check that
+ print *
+ print *,'IMPORTANT: it is your responsibility to make sure that in the input CUBIT (or similar) mesh'
+ print *,'that this code will read in SPECFEM3D format from files "nodes_coords_file" and "mesh_file"'
+ print *,'you have created layers of elements that constitute a layer of constant thickness aligned'
+ print *,'with the coordinate grid axes (X, Y and/or Z), so that this code can assign CPML flags to them.'
+ print *,'This code does NOT check that (because it cannot, in any easy way).'
+ print *,'The mesh inside these CPML layers does not need to be structured nor regular, any non-structured'
+ print *,'mesh is fine as long as it has flat PML inner and outer faces, parallel to the axes, and thus'
+ print *,'of a constant thickness.'
+ print *,'The thickness can be different between the X, Y and Z sides. But for X it must not vary,'
+ print *,'for Y it must not vary, and for Z it must not vary.'
+ print *,'If you do not know the exact thickness, or if you think the thickness may have tiny variations,'
+ print *,'for instance due to roundoff issues when exporting the mesh, you can use a slightly LARGER value'
+ print *,'in this code (say 2% to 5% more), and this code will fix that and will adjust it; never use'
+ print *,'a SMALLER value otherwise this code will miss some CPML elements.'
+ print *
+ print *,'Note: in a future release we will remove the constraint of having CPML layers aligned with the'
+ print *,'coordinate axes; we will allow for meshes that are titled by any constant angle in the horizontal plane.'
+ print *,'However this is not implemented yet.'
+ print *
+ print *,'1 = use a free surface at the top of the mesh (most classical option)'
+ print *,'2 = use a CPML absorbing layer at the top of the mesh (less classical option)'
+ print *,'3 = exit'
+ read(*,*) iflag
+ if(iflag /= 1 .and. iflag /= 2) stop 'exiting...'
+ if(iflag == 1) then
+ ALSO_ADD_ON_THE_TOP_SURFACE = .false.
+ else
+ ALSO_ADD_ON_THE_TOP_SURFACE = .true.
+ endif
+ print *
+
! open SPECFEM3D_Cartesian mesh file to read the points
open(unit=23,file='nodes_coords_file',status='old',action='read')
read(23,*) npoin
@@ -60,8 +89,40 @@
zmin = minval(z)
zmax = maxval(z)
-! ************* generate elements ******************
+ print *,'Xmin and Xmax of the mesh read = ',xmin,xmax
+ print *,'Ymin and Ymax of the mesh read = ',ymin,ymax
+ print *,'Zmin and Zmax of the mesh read = ',zmin,zmax
+ print *
+ print *,'What is the exact thickness of the PML layer that you want'
+ print *,'on the Xmin and Xmax faces of your mesh? (it needs to correspond exactly'
+ print *,'to the flat layers you created in your input CUBIT mesh, as mentioned in'
+ print *,'the comment printed above; if you think you have roundoff issues or very'
+ print *,'slightly varying thickness, give 2% or 5% more here, but never less'
+ read(*,*) THICKNESS_OF_X_PML
+ if(THICKNESS_OF_X_PML <= 0) stop 'negative thickness is not allowed; exiting...'
+ if(THICKNESS_OF_X_PML > 0.30*(xmax - xmin)) &
+ stop 'thickness of each CPML layer greater than 30% of the size of the mesh is not a good idea; exiting...'
+ print *
+
+ print *,'What is the exact thickness of the PML layer that you want'
+ print *,'on the Ymin and Ymax faces of your mesh?'
+ read(*,*) THICKNESS_OF_Y_PML
+ if(THICKNESS_OF_Y_PML <= 0) stop 'negative thickness is not allowed; exiting...'
+ if(THICKNESS_OF_Y_PML > 0.30*(ymax - ymin)) &
+ stop 'thickness of each CPML layer greater than 30% of the size of the mesh is not a good idea; exiting...'
+ print *
+
+ print *,'What is the exact thickness of the PML layer that you want'
+ print *,'on the Zmin and Zmax faces of your mesh?'
+ read(*,*) THICKNESS_OF_Z_PML
+ if(THICKNESS_OF_Z_PML <= 0) stop 'negative thickness is not allowed; exiting...'
+ if(THICKNESS_OF_Z_PML > 0.30*(zmax - zmin)) &
+ stop 'thickness of each CPML layer greater than 30% of the size of the mesh is not a good idea; exiting...'
+ print *
+
+! ************* read mesh elements and generate CPML flags *************
+
! open SPECFEM3D_Cartesian topology file to read the mesh elements
open(unit=23,file='mesh_file',status='old',action='read')
read(23,*) nspec
@@ -72,37 +133,14 @@
allocate(is_Y_CPML(nspec))
allocate(is_Z_CPML(nspec))
-! read local elements in this slice and output global DX elements
- do ispec=1,nspec
- read(23,*) ispec_read,i1,i2,i3,i4,i5,i6,i7,i8
- ibool(1,ispec_read) = i1
- ibool(2,ispec_read) = i2
- ibool(3,ispec_read) = i3
- ibool(4,ispec_read) = i4
- ibool(5,ispec_read) = i5
- ibool(6,ispec_read) = i6
- ibool(7,ispec_read) = i7
- ibool(8,ispec_read) = i8
- enddo
- close(23)
-
-! ************* generate CPML flags *************
-
is_X_CPML(:) = .false.
is_Y_CPML(:) = .false.
is_Z_CPML(:) = .false.
! loop on the whole mesh
- do ispec=1,nspec
+ do ispec_loop = 1,nspec
- i1 = ibool(1,ispec)
- i2 = ibool(2,ispec)
- i3 = ibool(3,ispec)
- i4 = ibool(4,ispec)
- i5 = ibool(5,ispec)
- i6 = ibool(6,ispec)
- i7 = ibool(7,ispec)
- i8 = ibool(8,ispec)
+ read(23,*) ispec,i1,i2,i3,i4,i5,i6,i7,i8
! Xmin CPML
limit = xmin + THICKNESS_OF_X_PML * SMALL_PERCENTAGE_TOLERANCE
@@ -138,6 +176,8 @@
enddo
+ close(23)
+
print *,'Total number of elements in the mesh = ',nspec
print *
print *,'Found ',count(is_X_CPML),' X_CPML elements'
@@ -164,10 +204,10 @@
! write the encoded CPML flag for each CPML element
do ispec=1,nspec
if(is_X_CPML(ispec) .and. is_Y_CPML(ispec) .and. is_Z_CPML(ispec)) then
- write(23,*) ispec,' 7'
+ write(23,*) ispec,' 7' !! DK DK mettre flag CPML_XYZ
else if(is_Y_CPML(ispec) .and. is_Z_CPML(ispec)) then
- write(23,*) ispec,' 6'
+ write(23,*) ispec,' 6' !! DK DK mettre flag CPML_YZ_ONLY
else if(is_X_CPML(ispec) .and. is_Z_CPML(ispec)) then
write(23,*) ispec,' 5'
@@ -188,5 +228,9 @@
enddo
close(23)
+ print *
+ print *,'CPML absorbing layer file "absorbing_cpml_file" has been successfully created'
+ print *
+
end program convert_mesh_to_CPML
More information about the CIG-COMMITS
mailing list