[cig-commits] r20764 - in seismo/2D/SPECFEM2D/trunk/src: meshfem2D specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Sat Sep 22 07:03:41 PDT 2012
Author: xie.zhinan
Date: 2012-09-22 07:03:40 -0700 (Sat, 22 Sep 2012)
New Revision: 20764
Modified:
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
add mpi support for using pml with external mesh
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90 2012-09-21 23:53:27 UTC (rev 20763)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90 2012-09-22 14:03:40 UTC (rev 20764)
@@ -374,6 +374,10 @@
! to store density and velocity model
integer, dimension(:), allocatable :: num_material
+ ! to store the position of pml element
+ integer, dimension(:), allocatable :: is_pml
+ integer :: nspec_cpml
+
! interface data
integer :: max_npoints_interface,number_of_interfaces,npoints_interface_bottom, &
npoints_interface_top
@@ -441,9 +445,13 @@
allocate(num_material(nelmnts))
num_material(:) = 0
+ allocate(is_pml(nelmnts))
+ is_pml(:) = 0
+
! assigns materials to mesh elements
if ( read_external_mesh ) then
call read_mat(materials_file, num_material)
+ call read_pml_element(CPML_element_file, is_pml, nspec_cpml)
else
call read_regions(nbregion,nb_materials,icodemat,cp,cs, &
rho_s,QKappa,Qmu,aniso3,aniso4,aniso5,aniso6,aniso7,aniso8, &
@@ -630,9 +638,9 @@
if ( any_abs ) then
call read_abs_surface(absorbing_surface_file, remove_min_to_start_at_zero)
-! rotate the elements that are located on the edges of the mesh if needed
-! otherwise the plane wave and Bielak conditions may not be applied correctly
- if(initialfield) call rotate_mesh_for_plane_wave(ngnod)
+ if(initialfield) then
+ call rotate_mesh_for_abs(ngnod)
+ endif
endif
else
@@ -968,7 +976,7 @@
endif
! *** generate the databases for the solver
- call save_databases(nspec,num_material, &
+ call save_databases(nspec,num_material, is_pml, &
my_interfaces,my_nb_interfaces, &
nnodes_tangential_curve,nodes_tangential_curve)
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/part_unstruct.F90 2012-09-21 23:53:27 UTC (rev 20763)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/part_unstruct.F90 2012-09-22 14:03:40 UTC (rev 20764)
@@ -201,7 +201,39 @@
end subroutine read_mat
+ !-----------------------------------------------
+ ! Read the position of pml element storing it in array 'is_pml'
+ !-----------------------------------------------
+ subroutine read_pml_element(filename, is_pml, nspec_cpml)
+ implicit none
+
+ character(len=256), intent(in) :: filename
+ integer, dimension(1:nelmnts), intent(out) :: is_pml
+ integer, intent(out) :: nspec_cpml
+! integer, dimension(:,:), allocatable :: local_pml
+
+ integer :: i,j,k,ier
+
+ open(unit=992, file=trim(filename), form='formatted' , status='old', action='read',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening file: ',trim(filename)
+ stop 'error read external mat file'
+ endif
+
+ read(992,*) nspec_cpml
+
+ do i = 1, nspec_cpml
+ read(992,*) j, k
+ is_pml(j) = k
+ enddo
+
+ close(992)
+
+ end subroutine read_pml_element
+
+
+
!-----------------------------------------------
! Read free surface.
! Edges from elastic elements are discarded.
@@ -330,13 +362,12 @@
end subroutine read_abs_surface
!-----------------------------------------------
- ! rotate_mesh_for_plane_wave.
- ! rotate mesh elements to make sure topological absorbing surface
- ! is aligned with physical absorbing surface, since this is the surface
- ! that we use to impose the plane wave and Bielak boundary conditions.
+ ! rotate_mesh_for_abs.
+ ! rotate mesh elements to make sure topological absorbing surface
+ ! is aligned with physical absorbing surface
!-----------------------------------------------
- subroutine rotate_mesh_for_plane_wave(ngnod)
+ subroutine rotate_mesh_for_abs(ngnod)
implicit none
@@ -344,8 +375,9 @@
integer :: i,j,ispec,i1,i2,inode,iswap
logical :: found_this_point
integer, dimension(:,:), allocatable :: ibool,ibool_rotated
-!! DK DK be careful here, "ibool" applies to the mesh corners (4 or 9 points) only,
+!! DK DK be careful here, "ibool" and "iglob" apply to the mesh corners (4 or 9 points) only,
!! DK DK not to the GLL points because there are no GLL points in the Gmsh mesh files
+! integer, dimension(:), allocatable :: ielem,iglob1,iglob2
integer :: index_rotation1,index_rotation2,index_rotation3,index_rotation4,&
index_rotation5,index_rotation6,index_rotation7,index_rotation8,index_edge
@@ -358,7 +390,7 @@
ibool(2,ispec) = elmnts((ispec-1)*ngnod+1)
ibool(3,ispec) = elmnts((ispec-1)*ngnod+2)
ibool(4,ispec) = elmnts((ispec-1)*ngnod+3)
- else if(ngnod == 9) then
+ elseif(ngnod == 9) then
ibool(1,ispec) = elmnts((ispec-1)*ngnod)
ibool(2,ispec) = elmnts((ispec-1)*ngnod+1)
ibool(3,ispec) = elmnts((ispec-1)*ngnod+2)
@@ -378,13 +410,13 @@
if(j==1)then
index_edge=3
ibool_rotated(:,:) = ibool(:,:)
- else if(j==2)then
+ elseif(j==2)then
index_edge=1
ibool(:,:) = ibool_rotated(:,:)
- else if(j==3)then
+ elseif(j==3)then
index_edge=4
ibool(:,:) = ibool_rotated(:,:)
- else if(j==4)then
+ elseif(j==4)then
index_edge=2
ibool(:,:) = ibool_rotated(:,:)
else
@@ -401,8 +433,8 @@
index_rotation6 = 4
index_rotation7 = 1
index_rotation8 = 4
- else if(index_edge==2)then
-! right edge
+ elseif(index_edge==2)then
+! right edge
index_rotation1 = 2
index_rotation2 = 3
index_rotation3 = 3
@@ -411,8 +443,8 @@
index_rotation6 = 4
index_rotation7 = 1
index_rotation8 = 2
- else if(index_edge==3)then
-! top edge
+ elseif(index_edge==3)then
+! top edge
index_rotation1 = 3
index_rotation2 = 4
index_rotation3 = 1
@@ -421,8 +453,8 @@
index_rotation6 = 2
index_rotation7 = 2
index_rotation8 = 3
- else if(index_edge==4)then
-! left edge
+ elseif(index_edge==4)then
+! left edge
index_rotation1 = 1
index_rotation2 = 4
index_rotation3 = 1
@@ -437,7 +469,7 @@
do i = 1,nelemabs
if(index_edge==abs_surface(5,i))then
- ispec = abs_surface(1,i) + 1 !!!! be careful: ispec from abs_surface(1,i) start at zero
+ ispec = abs_surface(1,i) + 1 !!!!be careful the ispec from abs_surface(1,i) begin from zero
found_this_point = .false.
do inode = 1,ngnod
if(ibool(inode,ispec) == abs_surface(3,i)) then
@@ -468,47 +500,46 @@
! test orientation
if(i1 == index_rotation1 .and. i2 == index_rotation2) then
-! print *,'orientation of element ',i,' is already good'
-
+! print *,'orientation of element ',i,' is already good'
else if (i1 == index_rotation3 .and. i2 == index_rotation4) then ! for this one, remember that we have swapped, thus 4 1 is 1 4
-! print *,'element ',i,' must be rotated 3 times'
+! print *,'element ',i,' must be rotated 3 times'
ibool_rotated(4,ispec) = ibool(1,ispec)
ibool_rotated(1,ispec) = ibool(2,ispec)
ibool_rotated(2,ispec) = ibool(3,ispec)
ibool_rotated(3,ispec) = ibool(4,ispec)
if(ngnod == 9) then
- ibool_rotated(8,ispec) = ibool(5,ispec)
- ibool_rotated(5,ispec) = ibool(6,ispec)
- ibool_rotated(6,ispec) = ibool(7,ispec)
- ibool_rotated(7,ispec) = ibool(8,ispec)
+ ibool_rotated(8,ispec) = ibool(5,ispec)
+ ibool_rotated(5,ispec) = ibool(6,ispec)
+ ibool_rotated(6,ispec) = ibool(7,ispec)
+ ibool_rotated(7,ispec) = ibool(8,ispec)
! 9th point is at the element center and thus never changes when we rotate an element
endif
else if (i1 == index_rotation5 .and. i2 == index_rotation6) then
-! print *,'element ',i,ispec,' must be rotated 2 times top'
+! print *,'element ',i,ispec,' must be rotated 2 times top'
ibool_rotated(3,ispec) = ibool(1,ispec)
ibool_rotated(4,ispec) = ibool(2,ispec)
ibool_rotated(1,ispec) = ibool(3,ispec)
ibool_rotated(2,ispec) = ibool(4,ispec)
if(ngnod == 9) then
- ibool_rotated(7,ispec) = ibool(5,ispec)
- ibool_rotated(8,ispec) = ibool(6,ispec)
- ibool_rotated(5,ispec) = ibool(7,ispec)
- ibool_rotated(6,ispec) = ibool(8,ispec)
+ ibool_rotated(7,ispec) = ibool(5,ispec)
+ ibool_rotated(8,ispec) = ibool(6,ispec)
+ ibool_rotated(5,ispec) = ibool(7,ispec)
+ ibool_rotated(6,ispec) = ibool(8,ispec)
! 9th point is at the element center and thus never changes when we rotate an element
endif
else if (i1 == index_rotation7 .and. i2 == index_rotation8) then
-! print *,'element ',i,' must be rotated 1 time'
+! print *,'element ',i,' must be rotated 1 time'
ibool_rotated(2,ispec) = ibool(1,ispec)
ibool_rotated(3,ispec) = ibool(2,ispec)
ibool_rotated(4,ispec) = ibool(3,ispec)
ibool_rotated(1,ispec) = ibool(4,ispec)
if(ngnod == 9) then
- ibool_rotated(6,ispec) = ibool(5,ispec)
- ibool_rotated(7,ispec) = ibool(6,ispec)
- ibool_rotated(8,ispec) = ibool(7,ispec)
- ibool_rotated(5,ispec) = ibool(8,ispec)
+ ibool_rotated(6,ispec) = ibool(5,ispec)
+ ibool_rotated(7,ispec) = ibool(6,ispec)
+ ibool_rotated(8,ispec) = ibool(7,ispec)
+ ibool_rotated(5,ispec) = ibool(8,ispec)
! 9th point is at the element center and thus never changes when we rotate an element
endif
@@ -528,7 +559,7 @@
elmnts((ispec-1)*ngnod+1) = ibool_rotated(2,ispec)
elmnts((ispec-1)*ngnod+2) = ibool_rotated(3,ispec)
elmnts((ispec-1)*ngnod+3) = ibool_rotated(4,ispec)
- else if(ngnod == 9) then
+ elseif(ngnod == 9) then
elmnts((ispec-1)*ngnod) = ibool_rotated(1,ispec)
elmnts((ispec-1)*ngnod+1) = ibool_rotated(2,ispec)
elmnts((ispec-1)*ngnod+2) = ibool_rotated(3,ispec)
@@ -543,7 +574,7 @@
endif
enddo
-end subroutine rotate_mesh_for_plane_wave
+end subroutine rotate_mesh_for_abs
!-----------------------------------------------
! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
@@ -810,7 +841,7 @@
! elastic element
is_acoustic_el = .false.
is_elastic_el = .true.
- else if ( phi_material(num_material(el+1)) >= 1.d0) then
+ elseif ( phi_material(num_material(el+1)) >= 1.d0) then
! acoustic element
is_acoustic_el = .true.
is_elastic_el = .false.
@@ -826,7 +857,7 @@
if ( phi_material(num_material(adjncy_g(el_adj)+1)) < TINYVAL) then
is_acoustic_el_adj = .false.
is_elastic_el_adj = .true.
- else if ( phi_material(num_material(adjncy_g(el_adj)+1)) >= 1.d0) then
+ elseif ( phi_material(num_material(adjncy_g(el_adj)+1)) >= 1.d0) then
is_acoustic_el_adj = .true.
is_elastic_el_adj = .false.
else
@@ -865,7 +896,7 @@
if ( phi_material(num_material(el+1)) < TINYVAL) then
is_acoustic_el = .false.
is_elastic_el = .true.
- else if ( phi_material(num_material(el+1)) >= 1.d0) then
+ elseif ( phi_material(num_material(el+1)) >= 1.d0) then
is_acoustic_el = .true.
is_elastic_el = .false.
else
@@ -876,7 +907,7 @@
if ( phi_material(num_material(adjncy_g(el_adj)+1)) < TINYVAL) then
is_acoustic_el_adj = .false.
is_elastic_el_adj = .true.
- else if ( phi_material(num_material(adjncy_g(el_adj)+1)) >= 1.d0) then
+ elseif ( phi_material(num_material(adjncy_g(el_adj)+1)) >= 1.d0) then
is_acoustic_el_adj = .true.
is_elastic_el_adj = .false.
else
@@ -959,7 +990,7 @@
! Write elements (their nodes) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
subroutine write_partition_database(IIN_database, iproc, nspec, &
- num_modele, ngnod, num_phase)
+ num_modele, num_pml, ngnod, num_phase)
implicit none
@@ -967,6 +998,7 @@
integer, intent(in) :: num_phase, iproc
integer, intent(inout) :: nspec
integer, dimension(:) :: num_modele
+ integer, dimension(:) :: num_pml
integer, intent(in) :: ngnod
integer :: i,j,k
@@ -989,7 +1021,7 @@
if (glob2loc_nodes_parts(k) == iproc) loc_nodes(j) = glob2loc_nodes(k)
enddo
enddo
- write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(i+1), (loc_nodes(k)+1, k=0,ngnod-1)
+ write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(i+1), (loc_nodes(k)+1, k=0,ngnod-1), num_pml(i+1)
endif
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90 2012-09-21 23:53:27 UTC (rev 20763)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/save_databases.f90 2012-09-22 14:03:40 UTC (rev 20764)
@@ -43,7 +43,7 @@
!========================================================================
- subroutine save_databases(nspec,num_material, &
+ subroutine save_databases(nspec,num_material,is_pml, &
my_interfaces,my_nb_interfaces, &
nnodes_tangential_curve,nodes_tangential_curve )
@@ -58,6 +58,7 @@
integer :: nspec
integer, dimension(nelmnts) :: num_material
+ integer, dimension(nelmnts) :: is_pml
integer, dimension(0:ninterfaces-1) :: my_interfaces
integer, dimension(0:ninterfaces-1) :: my_nb_interfaces
@@ -96,10 +97,11 @@
call write_glob2loc_nodes_database(15, iproc, npgeo, 1)
+! DK DK add support for using pml in mpi mode with external mesh
+! call write_partition_database(15, iproc, nspec, num_material, ngnod, 1)
+ call write_partition_database(15, iproc, nspec, num_material, is_pml, ngnod, 1)
- call write_partition_database(15, iproc, nspec, num_material, ngnod, 1)
-
write(15,*) 'npgeo nproc'
write(15,*) npgeo,nproc
@@ -123,7 +125,8 @@
if(read_external_mesh) then
write(15,*) 'CPML_element_file'
- write(15,"(a256)") trim(CPML_element_file)
+! write(15,"(a256)") trim(CPML_element_file)
+ write(15,*) trim(CPML_element_file)
endif
write(15,*) 'NELEM_PML_THICKNESS'
@@ -284,7 +287,9 @@
write(15,*) 'Arrays kmato and knods for each bloc:'
- call write_partition_database(15, iproc, nspec, num_material, ngnod, 2)
+! DK DK add support for using pml in mpi mode with external mesh
+! call write_partition_database(15, iproc, nspec, num_material, ngnod, 2)
+ call write_partition_database(15, iproc, nspec, num_material, is_pml, ngnod, 2)
if ( nproc /= 1 ) then
call write_interfaces_database(15, nproc, iproc, &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-09-21 23:53:27 UTC (rev 20763)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-09-22 14:03:40 UTC (rev 20764)
@@ -44,7 +44,7 @@
subroutine pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
- read_external_mesh,CPML_element_file)
+ read_external_mesh,region_CPML)
implicit none
@@ -70,10 +70,7 @@
!! DK DK for CPML_element_file
logical :: read_external_mesh
- character(len=256) :: CPML_element_file
- integer :: ier,ispec_CPML
- integer :: ispec_temp
- integer, dimension(:), allocatable :: region_CPML
+ integer, dimension(nspec) :: region_CPML
!!!detection of PML elements
@@ -178,60 +175,54 @@
is_PML(:) = .false.
which_PML_elem(:,:) = .false.
- open(unit=9999, file=trim(CPML_element_file), form='formatted' , status='old', action='read',iostat=ier)
- if( ier /= 0 ) then
- print*,'error opening file: ',trim(CPML_element_file)
- stop 'error read external CPML element file'
- endif
-
- read(9999,*)nspec_PML
- allocate(region_CPML(nspec_PML))
-
- do ispec_CPML=1,nspec_PML
- read(9999,*)ispec_temp,region_CPML(ispec_CPML)
- is_PML(ispec_temp)=.true.
- spec_to_PML(ispec_temp)=ispec_CPML
+ nspec_PML = 0
+ do ispec=1,nspec
+ if(region_CPML(ispec) .ne. 0)then
+ nspec_PML = nspec_PML + 1
+ is_PML(ispec)=.true.
+ spec_to_PML(ispec)=nspec_PML
+ endif
enddo
+
do ispec=1,nspec
if(is_PML(ispec))then
- ispec_CPML=spec_to_PML(ispec)
- if(region_CPML(ispec_CPML)==1)then
+ if(region_CPML(ispec)==1)then
which_PML_elem(ILEFT,ispec) = .true.
which_PML_elem(IRIGHT,ispec) = .false.
which_PML_elem(ITOP,ispec) = .false.
which_PML_elem(IBOTTOM,ispec) = .false.
- elseif(region_CPML(ispec_CPML)==2)then
+ elseif(region_CPML(ispec)==2)then
which_PML_elem(ILEFT,ispec) = .false.
which_PML_elem(IRIGHT,ispec) = .true.
which_PML_elem(ITOP,ispec) = .false.
which_PML_elem(IBOTTOM,ispec) = .false.
- elseif(region_CPML(ispec_CPML)==4)then
+ elseif(region_CPML(ispec)==4)then
which_PML_elem(ILEFT,ispec) = .false.
which_PML_elem(IRIGHT,ispec) = .false.
which_PML_elem(ITOP,ispec) = .true.
which_PML_elem(IBOTTOM,ispec) = .false.
- elseif(region_CPML(ispec_CPML)==5)then
+ elseif(region_CPML(ispec)==5)then
which_PML_elem(ILEFT,ispec) = .true.
which_PML_elem(IRIGHT,ispec) = .false.
which_PML_elem(ITOP,ispec) = .true.
which_PML_elem(IBOTTOM,ispec) = .false.
- elseif(region_CPML(ispec_CPML)==6)then
+ elseif(region_CPML(ispec)==6)then
which_PML_elem(ILEFT,ispec) = .false.
which_PML_elem(IRIGHT,ispec) = .true.
which_PML_elem(ITOP,ispec) = .true.
which_PML_elem(IBOTTOM,ispec) = .false.
- elseif(region_CPML(ispec_CPML)==8)then
+ elseif(region_CPML(ispec)==8)then
which_PML_elem(ILEFT,ispec) = .false.
which_PML_elem(IRIGHT,ispec) = .false.
which_PML_elem(ITOP,ispec) = .false.
which_PML_elem(IBOTTOM,ispec) = .true.
- elseif(region_CPML(ispec_CPML)==9)then
+ elseif(region_CPML(ispec)==9)then
which_PML_elem(ILEFT,ispec) = .true.
which_PML_elem(IRIGHT,ispec) = .false.
which_PML_elem(ITOP,ispec) = .false.
which_PML_elem(IBOTTOM,ispec) = .true.
- elseif(region_CPML(ispec_CPML)==10)then
+ elseif(region_CPML(ispec)==10)then
which_PML_elem(ILEFT,ispec) = .false.
which_PML_elem(IRIGHT,ispec) = .true.
which_PML_elem(ITOP,ispec) = .false.
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2012-09-21 23:53:27 UTC (rev 20763)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2012-09-22 14:03:40 UTC (rev 20764)
@@ -193,7 +193,8 @@
if(read_external_mesh)then
read(IIN,"(a80)") datlin
- read(IIN,"(a256)") CPML_element_file
+! read(IIN,"(a256)") CPML_element_file
+ read(IIN,*) CPML_element_file
endif
read(IIN,"(a80)") datlin
@@ -537,7 +538,7 @@
!
subroutine read_databases_mato(ipass,nspec,ngnod,kmato,knods, &
- perm,antecedent_list)
+ perm,antecedent_list,region_CPML)
! reads spectral macrobloc data
@@ -546,18 +547,20 @@
integer :: ipass,ngnod,nspec
integer, dimension(nspec) :: kmato
+ integer, dimension(nspec) :: region_CPML
integer, dimension(ngnod,nspec) :: knods
integer, dimension(nspec) :: perm,antecedent_list
! local parameters
- integer :: n,k,ispec,kmato_read
+ integer :: n,k,ispec,kmato_read,pml_read
integer, dimension(:), allocatable :: knods_read
character(len=80) :: datlin
! initializes
kmato(:) = 0
knods(:,:) = 0
+ region_CPML(:) = 0
! reads spectral macrobloc data
read(IIN,"(a80)") datlin
@@ -567,14 +570,16 @@
n = 0
do ispec = 1,nspec
! format: #element_id #material_id #node_id1 #node_id2 #...
- read(IIN,*) n,kmato_read,(knods_read(k), k=1,ngnod)
+ read(IIN,*) n,kmato_read,(knods_read(k), k=1,ngnod),pml_read
if(ipass == 1) then
! material association
kmato(n) = kmato_read
+ region_CPML(n) = pml_read
! element control node indices
knods(:,n)= knods_read(:)
else if(ipass == 2) then
kmato(perm(antecedent_list(n))) = kmato_read
+ region_CPML(perm(antecedent_list(n))) = pml_read
knods(:,perm(antecedent_list(n)))= knods_read(:)
else
call exit_MPI('error: maximum is 2 passes')
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-21 23:53:27 UTC (rev 20763)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-22 14:03:40 UTC (rev 20764)
@@ -999,7 +999,7 @@
!PML parameters
logical, dimension(:), allocatable :: is_PML
-
+ integer, dimension(:), allocatable :: region_CPML
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
@@ -1226,8 +1226,12 @@
allocate(antecedent_list(nspec))
allocate(perm(nspec))
endif
+
+ allocate(is_PML(nspec))
+! DK DK add support for using pml in mpi mode with external mesh
+ allocate(region_CPML(nspec))
call read_databases_mato(ipass,nspec,ngnod,kmato,knods, &
- perm,antecedent_list)
+ perm,antecedent_list,region_CPML)
!! DK DK Dec 2011: add a small crack (discontinuity) in the medium manually
if(ADD_A_SMALL_CRACK_IN_THE_MEDIUM) then
@@ -2850,19 +2854,18 @@
if( PML_BOUNDARY_CONDITIONS .and. anyabs_glob ) then
!PML code
- allocate(is_PML(nspec))
allocate(icorner_iglob(nglob))
allocate(which_PML_elem(4,nspec))
allocate(spec_to_PML(nspec))
is_PML(:) = .false.
which_PML_elem(:,:) = .false.
-
+! DK DK add support for using pml in mpi mode with external mesh
call pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
- read_external_mesh,CPML_element_file)
-
+ read_external_mesh,region_CPML)
+ deallocate(region_CPML)
deallocate(icorner_iglob)
if (nspec_PML==0) nspec_PML=1
More information about the CIG-COMMITS
mailing list