[cig-commits] [commit] devel: fix the problem that PMLs do not support the fluid-solid boundary inside CPML_Z_ONLY (aecbecc)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Dec 31 04:36:56 PST 2014
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/b259d542d122e904bd6a5f1b8519b8d46c7e0761...a7ad0d4491b67c622ec71c8e6297e93ff3bf4807
>---------------------------------------------------------------
commit aecbecc28586cab91b38b201d95a4942ef4ae02e
Author: Xie Zhinan <xiezhinan1984 at gmail.com>
Date: Wed Dec 31 19:02:18 2014 +0800
fix the problem that PMLs do not support the fluid-solid boundary inside CPML_Z_ONLY
>---------------------------------------------------------------
aecbecc28586cab91b38b201d95a4942ef4ae02e
list_for_Zhinan_for_SPECFEM2D_PMLs.txt | 75 ---------------------
src/specfem2D/compute_coupling_acoustic_el.f90 | 78 ++++++++++++----------
src/specfem2D/compute_coupling_viscoelastic_ac.f90 | 61 ++++++++---------
3 files changed, 71 insertions(+), 143 deletions(-)
diff --git a/list_for_Zhinan_for_SPECFEM2D_PMLs.txt b/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
index ef31622..4f9649a 100644
--- a/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
+++ b/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
@@ -1,80 +1,5 @@
For SPECFEM2D PMLs:
--------------------
-
-1/ current SPECFEM2D PMLs do not support the fluid-solid boundary not inside CPML_X_ONLY, see https://github.com/geodynamics/specfem2d/issues/81
-
-Date: Sat, 31 May 2014 23:50:43 +0200
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: è°¢å¿—å— <xiezhinan1984, Paul Cristini
-
-Dear Zhinan,
-
-There was also this one.
-If it is not too time consuming, it would be great to do it as well,
-this way all the minor issues detected in the last few months will
-finally be permanently fixed (this one (below) is the last one that I
-had noted; there a no other ones I think).
-
-Thanks a lot,
-Dimitri.
-
-Subject: Re: PML do not support the fluid-solid boundary not inside
-CPML_X_ONLY
-Date: Fri, 28 Feb 2014 18:51:18 +0100
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: è°¢å¿—å— <xiezhinan1984
-
-Dear Zhinan,
-
-It would be really great if you could do it if possible when you have
-time, because in the attached Par_file I had to rotate the model by 90
-degrees in order for the code to work (the model is supposed to be
-horizontal instead of vertical).
-
-Thanks again,
-Best wishes,
-
-Dimitri.
-
-On 26/02/2014 03:36, è°¢å¿—å— wrote:
-> Dear Dimitri,
->
-> We can to change to support CPML_Z_ONLY as well.
-> If you want, I can do the job.
-> Please tell the deadline that you think the code should be ready.
->
-> Thank you so much, Sir.
-> Best regards,
-> Zhinan
->
->
-> 2014-02-26 9:46 GMT+08:00 Dimitri Komatitsch:
->
-> Dear Zhinan,
->
-> When I try to put a fluid-solid boundary inside the PML in the
-> vertical direction, the code currently displays this:
->
-> PML do not support the fluid-solid boundary not inside CPML_X_ONLY.
->
-> How easy would it be to implement it in CPML_Z_ONLY as well?
->
-> (we could exclude the corners, i.e. change it to:
->
-> PML do not support the fluid-solid boundary inside CPML corners.
->
-> Could you tell me if you think this is easy to implement, or
-> difficult to implement?
->
-> Thank you very much,
-> Best regards,
->
-> Dimitri.
-
--------------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------------
2/ example "LuoYang_fluid_solid_kernel" of the 2D code seems to blow up when Stacey conditions are used:
diff --git a/src/specfem2D/compute_coupling_acoustic_el.f90 b/src/specfem2D/compute_coupling_acoustic_el.f90
index b0c00b8..43ac5c6 100644
--- a/src/specfem2D/compute_coupling_acoustic_el.f90
+++ b/src/specfem2D/compute_coupling_acoustic_el.f90
@@ -97,44 +97,52 @@
if(is_PML(ispec_elastic) .and. nspec_PML > 0) then
ispec_PML = spec_to_PML(ispec_elastic)
CPML_region_local = region_CPML(ispec_elastic)
+ kappa_x = K_x_store(i,j,ispec_PML)
+ kappa_z = K_z_store(i,j,ispec_PML)
+ d_x = d_x_store(i,j,ispec_PML)
+ d_z = d_z_store(i,j,ispec_PML)
+ alpha_x = alpha_x_store(i,j,ispec_PML)
+ alpha_z = alpha_z_store(i,j,ispec_PML)
+ beta_x = alpha_x + d_x / kappa_x
+ beta_z = alpha_z + d_z / kappa_z
+
if(CPML_region_local == CPML_X_ONLY)then
- kappa_x = K_x_store(i,j,ispec_PML)
- kappa_z = K_z_store(i,j,ispec_PML)
- d_x = d_x_store(i,j,ispec_PML)
- d_z = d_z_store(i,j,ispec_PML)
- alpha_x = alpha_x_store(i,j,ispec_PML)
- alpha_z = alpha_z_store(i,j,ispec_PML)
- beta_x = alpha_x + d_x / kappa_x
- beta_z = alpha_z + d_z / kappa_z
- call lik_parameter_computation(timeval,deltat,kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z,&
- CPML_region_local,13,A8,A9,A10,singularity_type_xz,bb_xz_1,bb_xz_2,&
- coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2)
- if(stage_time_scheme == 1) then
- rmemory_fsb_displ_elastic(1,1,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,1,i,j,inum) + &
- coef1_xz_1 * displ_elastic(1,iglob) + coef2_xz_1 * displ_elastic_old(1,iglob)
- rmemory_fsb_displ_elastic(1,3,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,3,i,j,inum) + &
- coef1_xz_1 * displ_elastic(3,iglob) + coef2_xz_1 * displ_elastic_old(3,iglob)
- endif
-
- if(stage_time_scheme == 6) then
- rmemory_fsb_displ_elastic_LDDRK(1,1,i,j,inum) = &
- alpha_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,1,i,j,inum) + &
- deltat * ( - bb_xz_1 * rmemory_fsb_displ_elastic(1,1,i,j,inum) + displ_elastic(1,iglob) )
- rmemory_fsb_displ_elastic(1,1,i,j,inum) = rmemory_fsb_displ_elastic(1,1,i,j,inum) + &
- beta_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,1,i,j,inum)
-
- rmemory_fsb_displ_elastic_LDDRK(1,3,i,j,inum) = &
- alpha_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,3,i,j,inum) + &
- deltat * ( - bb_xz_1 * rmemory_fsb_displ_elastic(1,3,i,j,inum) + displ_elastic(3,iglob) )
- rmemory_fsb_displ_elastic(1,3,i,j,inum) = rmemory_fsb_displ_elastic(1,3,i,j,inum) + &
- beta_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,3,i,j,inum)
- endif
-
- displ_x = A8 * displ_elastic(1,iglob) + A9 * rmemory_fsb_displ_elastic(1,1,i,j,inum)
- displ_z = A8 * displ_elastic(3,iglob) + A9 * rmemory_fsb_displ_elastic(1,3,i,j,inum)
+ call lik_parameter_computation(timeval,deltat,kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z,&
+ CPML_region_local,13,A8,A9,A10,singularity_type_xz,bb_xz_1,bb_xz_2,&
+ coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2)
+ elseif(CPML_region_local == CPML_Z_ONLY)then
+ call lik_parameter_computation(timeval,deltat,kappa_z,beta_z,alpha_z,kappa_x,beta_x,alpha_x,&
+ CPML_region_local,31,A8,A9,A10,singularity_type_xz,bb_xz_1,bb_xz_2,&
+ coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2)
+
else
- stop 'PML currently does not support a fluid-solid boundary located in a PML that is not CPML_X_ONLY'
+ stop 'PML do not support a fluid-solid boundary in corner PML region'
+ endif
+
+
+ if(stage_time_scheme == 1) then
+ rmemory_fsb_displ_elastic(1,1,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,1,i,j,inum) + &
+ coef1_xz_1 * displ_elastic(1,iglob) + coef2_xz_1 * displ_elastic_old(1,iglob)
+ rmemory_fsb_displ_elastic(1,3,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,3,i,j,inum) + &
+ coef1_xz_1 * displ_elastic(3,iglob) + coef2_xz_1 * displ_elastic_old(3,iglob)
+ endif
+
+ if(stage_time_scheme == 6) then
+ rmemory_fsb_displ_elastic_LDDRK(1,1,i,j,inum) = &
+ alpha_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,1,i,j,inum) + &
+ deltat * ( - bb_xz_1 * rmemory_fsb_displ_elastic(1,1,i,j,inum) + displ_elastic(1,iglob) )
+ rmemory_fsb_displ_elastic(1,1,i,j,inum) = rmemory_fsb_displ_elastic(1,1,i,j,inum) + &
+ beta_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,1,i,j,inum)
+
+ rmemory_fsb_displ_elastic_LDDRK(1,3,i,j,inum) = &
+ alpha_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,3,i,j,inum) + &
+ deltat * ( - bb_xz_1 * rmemory_fsb_displ_elastic(1,3,i,j,inum) + displ_elastic(3,iglob) )
+ rmemory_fsb_displ_elastic(1,3,i,j,inum) = rmemory_fsb_displ_elastic(1,3,i,j,inum) + &
+ beta_LDDRK(i_stage) * rmemory_fsb_displ_elastic_LDDRK(1,3,i,j,inum)
endif
+
+ displ_x = A8 * displ_elastic(1,iglob) + A9 * rmemory_fsb_displ_elastic(1,1,i,j,inum)
+ displ_z = A8 * displ_elastic(3,iglob) + A9 * rmemory_fsb_displ_elastic(1,3,i,j,inum)
else
displ_x = displ_elastic(1,iglob)
displ_z = displ_elastic(3,iglob)
diff --git a/src/specfem2D/compute_coupling_viscoelastic_ac.f90 b/src/specfem2D/compute_coupling_viscoelastic_ac.f90
index 0c7e90b..9fc5c59 100644
--- a/src/specfem2D/compute_coupling_viscoelastic_ac.f90
+++ b/src/specfem2D/compute_coupling_viscoelastic_ac.f90
@@ -95,40 +95,35 @@
if(is_PML(ispec_acoustic) .and. nspec_PML > 0) then
ispec_PML = spec_to_PML(ispec_acoustic)
CPML_region_local = region_CPML(ispec_acoustic)
- if(CPML_region_local == CPML_X_ONLY)then
- kappa_x = K_x_store(i,j,ispec_PML)
- kappa_z = K_z_store(i,j,ispec_PML)
- d_x = d_x_store(i,j,ispec_PML)
- d_z = d_z_store(i,j,ispec_PML)
- alpha_x = alpha_x_store(i,j,ispec_PML)
- alpha_z = alpha_z_store(i,j,ispec_PML)
- beta_x = alpha_x + d_x / kappa_x
- beta_z = alpha_z + d_z / kappa_z
- call l_parameter_computation(timeval,deltat,kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z, &
- CPML_region_local,A0,A1,A2,A3,A4,singularity_type,&
- bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2)
-
- if(stage_time_scheme == 1) then
- rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) = &
- coef0_1 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + &
- coef1_1 * potential_acoustic(iglob) + coef2_1 * potential_acoustic_old(iglob)
- endif
-
- if(stage_time_scheme == 6) then
- rmemory_sfb_potential_ddot_acoustic_LDDRK(1,i,j,inum) = &
- alpha_LDDRK(i_stage) * rmemory_sfb_potential_ddot_acoustic_LDDRK(1,i,j,inum) + &
- deltat * (-bb_1 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + potential_acoustic(iglob))
- rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) = rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + &
- beta_LDDRK(i_stage) * rmemory_sfb_potential_ddot_acoustic_LDDRK(1,i,j,inum)
-
- endif
-
- pressure = - (A0 * potential_dot_dot_acoustic(iglob) + A1 * potential_dot_acoustic(iglob) + &
- A2 * potential_acoustic(iglob) + A3 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum))
-
- else
- stop 'PML do not support the fluid-solid boundary not inside CPML_X_ONLY'
+ kappa_x = K_x_store(i,j,ispec_PML)
+ kappa_z = K_z_store(i,j,ispec_PML)
+ d_x = d_x_store(i,j,ispec_PML)
+ d_z = d_z_store(i,j,ispec_PML)
+ alpha_x = alpha_x_store(i,j,ispec_PML)
+ alpha_z = alpha_z_store(i,j,ispec_PML)
+ beta_x = alpha_x + d_x / kappa_x
+ beta_z = alpha_z + d_z / kappa_z
+ call l_parameter_computation(timeval,deltat,kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z, &
+ CPML_region_local,A0,A1,A2,A3,A4,singularity_type,&
+ bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2)
+
+ if(stage_time_scheme == 1) then
+ rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) = &
+ coef0_1 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + &
+ coef1_1 * potential_acoustic(iglob) + coef2_1 * potential_acoustic_old(iglob)
endif
+
+ if(stage_time_scheme == 6) then
+ rmemory_sfb_potential_ddot_acoustic_LDDRK(1,i,j,inum) = &
+ alpha_LDDRK(i_stage) * rmemory_sfb_potential_ddot_acoustic_LDDRK(1,i,j,inum) + &
+ deltat * (-bb_1 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + potential_acoustic(iglob))
+ rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) = rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + &
+ beta_LDDRK(i_stage) * rmemory_sfb_potential_ddot_acoustic_LDDRK(1,i,j,inum)
+
+ endif
+
+ pressure = - (A0 * potential_dot_dot_acoustic(iglob) + A1 * potential_dot_acoustic(iglob) + &
+ A2 * potential_acoustic(iglob) + A3 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum))
else
pressure = - potential_dot_dot_acoustic(iglob)
endif
More information about the CIG-COMMITS
mailing list