[cig-commits] [commit] devel: fix the current way for excluding corners to make sure there is no contradiction on the normal when Stacey ABC is used, which did not work, since for different IEDGE array "codeabs(IEDGE, ispecabs)" is independent of each other. (29a21a0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Dec 31 02:42:05 PST 2014
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/ed7588a2e02bfe450c867c17c4ee9701cfc69f45...b259d542d122e904bd6a5f1b8519b8d46c7e0761
>---------------------------------------------------------------
commit 29a21a0e12c7ea359fba4fe922f99fa13206a8df
Author: Xie Zhinan <xiezhinan1984 at gmail.com>
Date: Wed Dec 31 13:38:44 2014 +0800
fix the current way for excluding corners to make sure there is no contradiction on the normal when Stacey ABC is used,
which did not work, since for different IEDGE array "codeabs(IEDGE,ispecabs)" is independent of each other.
>---------------------------------------------------------------
29a21a0e12c7ea359fba4fe922f99fa13206a8df
.../plot_compare_to_analytical_solution.gnu | 16 +-
list_for_Zhinan_for_SPECFEM2D_PMLs.txt | 182 +++++----------------
src/specfem2D/compute_forces_acoustic.f90 | 10 +-
src/specfem2D/compute_forces_poro_fluid.f90 | 10 +-
src/specfem2D/compute_forces_poro_solid.f90 | 10 +-
src/specfem2D/compute_forces_viscoelastic.F90 | 6 +-
src/specfem2D/invert_mass_matrix.F90 | 14 +-
src/specfem2D/read_databases.F90 | 50 +++++-
src/specfem2D/specfem2D.F90 | 6 +
src/specfem2D/specfem2D_par.f90 | 4 +
10 files changed, 136 insertions(+), 172 deletions(-)
diff --git a/EXAMPLES/semi_infinite_homo/plot_compare_to_analytical_solution.gnu b/EXAMPLES/semi_infinite_homo/plot_compare_to_analytical_solution.gnu
index 0927a4b..f22d92a 100644
--- a/EXAMPLES/semi_infinite_homo/plot_compare_to_analytical_solution.gnu
+++ b/EXAMPLES/semi_infinite_homo/plot_compare_to_analytical_solution.gnu
@@ -6,39 +6,39 @@ set ylabel "Amplitude of displacement component (m)"
set xrange [-0.15:0.7]
-plot "OUTPUT_FILES/S0001.AA.BXX.semd" t 'Numerical Ux' w l lc 1
+plot "OUTPUT_FILES/AA.S0001.BXX.semd" t 'Numerical Ux' w l lc 1
pause -1 "Hit any key..."
plot "Ux_file.dat" t 'Analytic Ux' w l lc 3
pause -1 "Hit any key..."
-plot "OUTPUT_FILES/S0001.AA.BXX.semd" t 'Numerical Ux' w l lc 1, "Ux_file.dat" t 'Analytic Ux' w l lc 3
+plot "OUTPUT_FILES/AA.S0001.BXX.semd" t 'Numerical Ux' w l lc 1, "Ux_file.dat" t 'Analytic Ux' w l lc 3
pause -1 "Hit any key..."
-plot "OUTPUT_FILES/S0001.AA.BXZ.semd" t 'Numerical Uz' w l lc 1
+plot "OUTPUT_FILES/AA.S0001.BXZ.semd" t 'Numerical Uz' w l lc 1
pause -1 "Hit any key..."
plot "Uz_file.dat" t 'Analytic Uz' w l lc 3
pause -1 "Hit any key..."
-plot "OUTPUT_FILES/S0001.AA.BXZ.semd" t 'Numerical Uz' w l lc 1, "Uz_file.dat" t 'Analytic Uz' w l lc 3
+plot "OUTPUT_FILES/AA.S0001.BXZ.semd" t 'Numerical Uz' w l lc 1, "Uz_file.dat" t 'Analytic Uz' w l lc 3
pause -1 "Hit any key..."
#plot "Ux_file.dat" t 'Analytic Ux' w l lc 3
#pause -1 "Hit any key..."
-#plot "OUTPUT_FILES/S0001.AA.BXX.semd" t 'Numerical Ux' w l lc 1
+#plot "OUTPUT_FILES/AA.S0001.BXX.semd" t 'Numerical Ux' w l lc 1
#pause -1 "Hit any key..."
-#plot "OUTPUT_FILES/S0001.AA.BXZ.semd" t 'Numerical Uz' w l lc 1
+#plot "OUTPUT_FILES/AA.S0001.BXZ.semd" t 'Numerical Uz' w l lc 1
#pause -1 "Hit any key..."
#plot "source.txt" t 'Numerical Uz' w l lc 1
#pause -1 "Hit any key..."
-#plot "OUTPUT_FILES/S0001.AA.BXX.semd" t 'Numerical Ux' w l lc 1, "S0001.AA.BXX.rk.semd" t 'RKUx' w l lc 3
+#plot "OUTPUT_FILES/AA.S0001.BXX.semd" t 'Numerical Ux' w l lc 1, "AA.S0001.BXX.rk.semd" t 'RKUx' w l lc 3
#pause -1 "Hit any key..."
#
-#plot "OUTPUT_FILES/S0001.AA.BXZ.semd" t 'Numerical Uz' w l lc 1, "S0001.AA.BXZ.rk.semd" t 'RKUw' w l lc 3
+#plot "OUTPUT_FILES/AA.S0001.BXZ.semd" t 'Numerical Uz' w l lc 1, "AA.S0001.BXZ.rk.semd" t 'RKUw' w l lc 3
#pause -1 "Hit any key..."
diff --git a/list_for_Zhinan_for_SPECFEM2D_PMLs.txt b/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
index 2db6d54..ef31622 100644
--- a/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
+++ b/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
@@ -95,26 +95,6 @@ Best regards,
Dimitri.
-On 04/09/2014 04:45, xiezhinan wrote:
-> Dear Dimitri and Ying,
->
-> An initial tests show the problem is that:
->
-> Even for forward simulation with a purely fluid model, when turning the
-> bottom and top Stacey BC on, the simulation quickly blow up.
-> But for tests with right and left Stacey BC only, the simulation is fine.
->
-> While for forward simulation with a purely solid model,everything is fine.
->
-> Thus, I will first focus on solving the following problem:
-> Fixing Stacey BC for a purely fluid model with Stacey BC imposed on its
-> four out-edge boundaries.
->
-> Thank you so much.
->
-> Best regards,
-> Zhinan
->
>
> At 2014-09-04 05:56:50, "Dimitri Komatitsch" wrote:
>>
@@ -396,22 +376,15 @@ Zhinan
Hi Zhinan,
- I agree with you that the process of selecting the PML
- parameters
- should be as automatic as possible, i.e. if possible the
- code itself
- should decide and the user should not have to enter all the
- damping
+ I agree with you that the process of selecting the PML parameters
+ should be as automatic as possible, i.e. if possible the code itself
+ should decide and the user should not have to enter all the damping
parameters.
- However in the automatic decisions made by the code I think
- it makes
- sense to distinguish the thickness and the damping
- parameters along
- X from those along Z, since it is very common to have very
- elongated
- models (for instance in ocean acoustics, but also in oil
- industry
+ However in the automatic decisions made by the code I think it makes
+ sense to distinguish the thickness and the damping parameters along
+ X from those along Z, since it is very common to have very elongated
+ models (for instance in ocean acoustics, but also in oil industry
large-offset models).
Thanks a lot,
@@ -424,20 +397,12 @@ Zhinan
Hi Dimitri,
Thank you so much.
- After I get a more clean solution to Paul's problem, I
- will
- commit the
- change but set the streching function in different
- direction
- equal in
- default and leave a comment.
- Due to the fact that there is no clear guideline on
- choosing the
- parameters inside streching function, maybe it is
- useless to
- leave the
- general user a flexible choice to define the streching
- function.
+ After I get a more clean solution to Paul's problem, I will commit the
+ change but set the streching function in different direction
+ equal in default and leave a comment.
+ Due to the fact that there is no clear guideline on choosing the
+ parameters inside streching function, maybe it is useless to
+ leave the general user a flexible choice to define the streching function.
If not, please correct me, Sir.
Thanks again.
Best regards,
@@ -447,18 +412,12 @@ Zhinan
Hi Zhinan,
- Thank you so much. Could you please commit the
- changes to GIT?
-
- (if you think these changes are only useful for
- Paul but not in
- general, then please commit them as comments, i.e.
- leave
- the current
- values but put the new lines (the modified lines) as
- comments in the
- routines, and then please create a Pull Request on
- GitHub.
+ Thank you so much. Could you please commit the changes to GIT?
+
+ (if you think these changes are only useful for Paul but not in
+ general, then please commit them as comments, i.e. leave the current
+ values but put the new lines (the modified lines) as comments in the
+ routines, and then please create a Pull Request on GitHub.
(otherwise the changes will quickly be lost)
Thank you very much,
@@ -471,47 +430,21 @@ Zhinan
Hi Paul,
- I confirmed this problem comes from the grazing
- incident of wave.
-
- As suggested in our solid PML paper, the idea is
- originally from the
- paper of Zhang, we need to carefully design
- the stretching
- function to
- achieve an excellent absorption. I.E. we need
- to carefully
- choose both
- d, k, and alpha.
-
- I have tuned the parameter inside pml_init.F90 to
- achieve a good
- result
- for the example you send to me. Both the
- result and
- pml_init.F90
- have
- been attached.
-
- Currently the solution is not perfect, since I
- still
- have small
- reflection from the side PML due to the use of
- varying
- K, which is
- bigger than 1.
- But I think this problem can be solved by
- admitting
- different
- choice of
- d_max, k_max,, and alpha_max, along difference
- space
- coordinate,
- which
- is not allowed in our current implementation.
- I will
- keep you
- confirmed.
+ I confirmed this problem comes from the grazing incident of wave.
+
+ As suggested in our solid PML paper, the idea is originally from the
+ paper of Zhang, we need to carefully design the stretching function to
+ achieve an excellent absorption. I.E. we need to carefully
+ choose both d, k, and alpha.
+
+ I have tuned the parameter inside pml_init.F90 to achieve a good
+ result for the example you send to me. Both the result and pml_init.F90 have been attached.
+
+ Currently the solution is not perfect, since I still have small
+ reflection from the side PML due to the use of varying K, which is bigger than 1.
+ But I think this problem can be solved by admitting different choice of
+ d_max, k_max,, and alpha_max, along difference space coordinate, which
+ is not allowed in our current implementation. I will keep you confirmed.
Thanks again,
Best regards,
@@ -525,17 +458,12 @@ Zhinan
I am really sorry that I missed this email.
- The reflection come from the bottom PML,
- but seems
- work
- fine for
+ The reflection come from the bottom PML, but seems work fine for
side PML and the corner PML.
- Thus the problem may come from the wrong
- definition of the PML
+ Thus the problem may come from the wrong definition of the PML
element in the bottom
- I will have a test today and keep you
- confirmed.
+ I will have a test today and keep you confirmed.
Best regards,
Zhinan
@@ -545,34 +473,14 @@ Zhinan
Dear Zhinan,
- I am in the process of adding a new
- example in the EXAMPLES
- directory corresponding of our Jasa
- Express
- Letter on
- underwater
- acoustics and I am facing a problem
- with CPML
- in fluid
- domains.
- The configuration is a fluid layer
- overlying a
- sediment
- which
- can be either fluid or elastic but
- with the same
- velocity for P
- waves. As you can see in the attached
- file the
- CMPL is
- working
- when the sediment is elastic but when
- it is
- fluid there
- is a
- spurious reflection.
- Do you have any ideas of what's going
- on ?
+ I am in the process of adding a new example in the EXAMPLES
+ directory corresponding of our Jasa Express Letter on underwater
+ acoustics and I am facing a problem with CPML in fluid domains.
+ The configuration is a fluid layer overlying a sediment which
+ can be either fluid or elastic but with the same velocity for P
+ waves. As you can see in the attached file the CMPL is working
+ when the sediment is elastic but when it is fluid there is a spurious reflection.
+ Do you have any ideas of what's going on ?
Thanks
diff --git a/src/specfem2D/compute_forces_acoustic.f90 b/src/specfem2D/compute_forces_acoustic.f90
index f8bc468..f6fc9c6 100644
--- a/src/specfem2D/compute_forces_acoustic.f90
+++ b/src/specfem2D/compute_forces_acoustic.f90
@@ -49,7 +49,7 @@
use specfem_par, only: nglob,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs,acoustic, &
- codeabs,stage_time_scheme,i_stage, &
+ codeabs,codeabs_corner,stage_time_scheme,i_stage, &
density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
vpext,rhoext, &
hprime_xx,hprimewgll_xx, &
@@ -579,8 +579,8 @@
ibegin = ibegin_edge1(ispecabs)
iend = iend_edge1(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(1,ispecabs)) ibegin = 2
+ if(codeabs_corner(2,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
iglob = ibool(i,j,ispec)
! external velocity model
@@ -622,8 +622,8 @@
ibegin = ibegin_edge3(ispecabs)
iend = iend_edge3(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(3,ispecabs)) ibegin = 2
+ if(codeabs_corner(4,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
iglob = ibool(i,j,ispec)
! external velocity model
diff --git a/src/specfem2D/compute_forces_poro_fluid.f90 b/src/specfem2D/compute_forces_poro_fluid.f90
index 96d6fed..86f4862 100644
--- a/src/specfem2D/compute_forces_poro_fluid.f90
+++ b/src/specfem2D/compute_forces_poro_fluid.f90
@@ -49,7 +49,7 @@
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
- ibool,kmato,numabs,poroelastic,codeabs, &
+ ibool,kmato,numabs,poroelastic,codeabs,codeabs_corner, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
displs_poroelastic_old,b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
@@ -743,8 +743,8 @@
iend = iend_edge1_poro(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(1,ispecabs)) ibegin = 2
+ if(codeabs_corner(2,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
@@ -803,8 +803,8 @@
iend = iend_edge3_poro(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(3,ispecabs)) ibegin = 2
+ if(codeabs_corner(4,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
diff --git a/src/specfem2D/compute_forces_poro_solid.f90 b/src/specfem2D/compute_forces_poro_solid.f90
index 44686d5..cb43cf2 100644
--- a/src/specfem2D/compute_forces_poro_solid.f90
+++ b/src/specfem2D/compute_forces_poro_solid.f90
@@ -49,7 +49,7 @@
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
- ibool,kmato,numabs,poroelastic,codeabs, &
+ ibool,kmato,numabs,poroelastic,codeabs,codeabs_corner, &
accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displs_poroelastic_old,&
displw_poroelastic,b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
@@ -764,8 +764,8 @@
iend = iend_edge1_poro(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(1,ispecabs)) ibegin = 2
+ if(codeabs_corner(2,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
@@ -824,8 +824,8 @@
iend = iend_edge3_poro(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(3,ispecabs)) ibegin = 2
+ if(codeabs_corner(4,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
diff --git a/src/specfem2D/compute_forces_viscoelastic.F90 b/src/specfem2D/compute_forces_viscoelastic.F90
index ff668df..fe6ed15 100644
--- a/src/specfem2D/compute_forces_viscoelastic.F90
+++ b/src/specfem2D/compute_forces_viscoelastic.F90
@@ -52,7 +52,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
- ibool,kmato,numabs,elastic,codeabs, &
+ ibool,kmato,numabs,elastic,codeabs,codeabs_corner, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext,&
source_time_function,sourcearray,adj_sourcearrays,anisotropic,anisotropy, &
@@ -1355,7 +1355,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic
! exclude corners to make sure there is no contradiction on the normal
! for Stacey absorbing conditions but not for incident plane waves;
! thus subtract nothing i.e. zero in that case
- if((codeabs(IEDGE4,ispecabs) .and. i == 1) .or. (codeabs(IEDGE2,ispecabs) .and. i == NGLLX)) then
+ if((codeabs_corner(1,ispecabs) .and. i == 1) .or. (codeabs_corner(2,ispecabs) .and. i == NGLLX)) then
tx = 0._CUSTOM_REAL; ty = 0._CUSTOM_REAL; tz = 0._CUSTOM_REAL
endif
@@ -1465,7 +1465,7 @@ subroutine compute_forces_viscoelastic(accel_elastic,veloc_elastic,displ_elastic
! exclude corners to make sure there is no contradiction on the normal
! for Stacey absorbing conditions but not for incident plane waves;
! thus subtract nothing i.e. zero in that case
- if((codeabs(IEDGE4,ispecabs) .and. i == 1) .or. (codeabs(IEDGE2,ispecabs) .and. i == NGLLX)) then
+ if((codeabs_corner(3,ispecabs) .and. i == 1) .or. (codeabs_corner(4,ispecabs) .and. i == NGLLX)) then
tx = 0._CUSTOM_REAL; ty = 0._CUSTOM_REAL; tz = 0._CUSTOM_REAL
endif
diff --git a/src/specfem2D/invert_mass_matrix.F90 b/src/specfem2D/invert_mass_matrix.F90
index 58a829f..96a0d84 100644
--- a/src/specfem2D/invert_mass_matrix.F90
+++ b/src/specfem2D/invert_mass_matrix.F90
@@ -58,7 +58,7 @@
assign_external_model,numat, &
density,poroelastcoef,porosity,tortuosity, &
vpext,rhoext, &
- anyabs,numabs,deltat,codeabs,&
+ anyabs,numabs,deltat,codeabs,codeabs_corner,&
ibegin_edge1,iend_edge1,ibegin_edge3,iend_edge3, &
ibegin_edge4,iend_edge4,ibegin_edge2,iend_edge2, &
rmass_inverse_elastic_three,&
@@ -513,7 +513,7 @@
! exclude corners to make sure there is no contradiction on the normal
! for Stacey absorbing conditions but not for incident plane waves;
! thus subtract nothing i.e. zero in that case
- if((codeabs(IEDGE4,ispecabs) .and. i == 1) .or. (codeabs(IEDGE2,ispecabs) .and. i == NGLLX)) then
+ if((codeabs_corner(1,ispecabs) .and. i == 1) .or. (codeabs_corner(2,ispecabs) .and. i == NGLLX)) then
tx = 0
ty = 0
tz = 0
@@ -574,7 +574,7 @@
! exclude corners to make sure there is no contradiction on the normal
! for Stacey absorbing conditions but not for incident plane waves;
! thus subtract nothing i.e. zero in that case
- if((codeabs(IEDGE4,ispecabs) .and. i == 1) .or. (codeabs(IEDGE2,ispecabs) .and. i == NGLLX)) then
+ if((codeabs_corner(3,ispecabs) .and. i == 1) .or. (codeabs_corner(4,ispecabs) .and. i == NGLLX)) then
tx = 0
ty = 0
tz = 0
@@ -660,8 +660,8 @@
ibegin = ibegin_edge1(ispecabs)
iend = iend_edge1(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(1,ispecabs)) ibegin = 2
+ if(codeabs_corner(2,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
iglob = ibool(i,j,ispec)
! external velocity model
@@ -686,8 +686,8 @@
ibegin = ibegin_edge3(ispecabs)
iend = iend_edge3(ispecabs)
! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(IEDGE4,ispecabs)) ibegin = 2
- if(codeabs(IEDGE2,ispecabs)) iend = NGLLX-1
+ if(codeabs_corner(3,ispecabs)) ibegin = 2
+ if(codeabs_corner(4,ispecabs)) iend = NGLLX-1
do i = ibegin,iend
iglob = ibool(i,j,ispec)
! external velocity model
diff --git a/src/specfem2D/read_databases.F90 b/src/specfem2D/read_databases.F90
index f4336e9..643daec 100644
--- a/src/specfem2D/read_databases.F90
+++ b/src/specfem2D/read_databases.F90
@@ -598,7 +598,7 @@
use specfem_par, only : myrank,nelemabs,nspec,anyabs, &
ibegin_edge1,iend_edge1,ibegin_edge2,iend_edge2, &
ibegin_edge3,iend_edge3,ibegin_edge4,iend_edge4, &
- numabs,codeabs,typeabs, &
+ numabs,codeabs,codeabs_corner,typeabs, &
nspec_left,nspec_right,nspec_bottom,nspec_top, &
ib_right,ib_left,ib_bottom,ib_top,PML_BOUNDARY_CONDITIONS
@@ -606,7 +606,7 @@
include "constants.h"
! local parameters
- integer :: inum,numabsread,typeabsread
+ integer :: inum,inum_duplicate,numabsread,typeabsread
logical :: codeabsread(4)
character(len=80) :: datlin
@@ -618,6 +618,7 @@
! initializes
codeabs(:,:) = .false.
+ codeabs_corner(:,:) = .false.
typeabs(:) = 0
ibegin_edge1(:) = 0
@@ -679,6 +680,51 @@
enddo
+! detection of the corner element
+ do inum = 1,nelemabs
+
+ if (codeabs(IEDGE1,inum)) then
+ do inum_duplicate = 1,nelemabs
+ if (inum == inum_duplicate) then
+ ! left for blank, since no operation is needed.
+ else
+ if (numabs(inum) == numabs(inum_duplicate)) then
+ if (codeabs(IEDGE4,inum_duplicate)) then
+ codeabs_corner(1,inum) = .true.
+ endif
+
+ if (codeabs(IEDGE2,inum_duplicate)) then
+ codeabs_corner(2,inum) = .true.
+ endif
+
+ endif
+ endif
+ enddo
+ endif
+
+ if (codeabs(IEDGE3,inum)) then
+ do inum_duplicate = 1,nelemabs
+ if (inum == inum_duplicate) then
+ ! left for blank, since no operation is needed.
+ else
+ if (numabs(inum) == numabs(inum_duplicate)) then
+ if (codeabs(IEDGE4,inum_duplicate)) then
+ codeabs_corner(3,inum) = .true.
+ endif
+
+ if (codeabs(IEDGE2,inum_duplicate)) then
+ codeabs_corner(4,inum) = .true.
+ endif
+
+ endif
+ endif
+ enddo
+ endif
+
+ enddo
+
+! detection of the corner element
+
! boundary element numbering
do inum = 1,nelemabs
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index 3e3ce09..af406e4 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -802,6 +802,12 @@
allocate(numabs(nelemabs))
allocate(codeabs(4,nelemabs))
+
+!---codeabs_corner(1,nelemabs) denotes whether element is on bottom-left corner of absorbing boundary or not
+!---codeabs_corner(2,nelemabs) denotes whether element is on bottom-right corner of absorbing boundary or not
+!---codeabs_corner(3,nelemabs) denotes whether element is on top-left corner of absorbing boundary or not
+!---codeabs_corner(4,nelemabs) denotes whether element is on top-right corner of absorbing boundary or not
+ allocate(codeabs_corner(4,nelemabs))
allocate(typeabs(nelemabs))
allocate(ibegin_edge1(nelemabs))
diff --git a/src/specfem2D/specfem2D_par.f90 b/src/specfem2D/specfem2D_par.f90
index e6179ae..7785bb1 100644
--- a/src/specfem2D/specfem2D_par.f90
+++ b/src/specfem2D/specfem2D_par.f90
@@ -305,6 +305,10 @@ module specfem_par
logical, dimension(:,:), allocatable :: codeabs
integer, dimension(:), allocatable :: typeabs
+! for detection of corner element on absorbing boundary
+ logical, dimension(:,:), allocatable :: codeabs_corner
+
+
! for attenuation
integer :: N_SLS
double precision, dimension(:), allocatable :: QKappa_attenuation
More information about the CIG-COMMITS
mailing list