[cig-commits] [commit] devel: remove the staff in the list_for_Zhinan_for_SPECFEM3D_PMLs.txt concerning the implementation of elastic PML for viscoelastic wave simulation (3b04b68)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Jan 27 03:36:44 PST 2015
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/367727d9b25a3dc40b9c10fee9dc3275eff7b717...3b04b68452f55959366a6a10017252704d6e711d
>---------------------------------------------------------------
commit 3b04b68452f55959366a6a10017252704d6e711d
Author: Xie Zhinan <xiezhinan1984 at gmail.com>
Date: Tue Jan 27 06:31:06 2015 +0100
remove the staff in the list_for_Zhinan_for_SPECFEM3D_PMLs.txt
concerning the implementation of elastic PML for viscoelastic wave simulation
>---------------------------------------------------------------
3b04b68452f55959366a6a10017252704d6e711d
list_for_Zhinan_for_SPECFEM3D_PMLs.txt | 268 ---------------------------------
1 file changed, 268 deletions(-)
diff --git a/list_for_Zhinan_for_SPECFEM3D_PMLs.txt b/list_for_Zhinan_for_SPECFEM3D_PMLs.txt
index a470553..71c9465 100644
--- a/list_for_Zhinan_for_SPECFEM3D_PMLs.txt
+++ b/list_for_Zhinan_for_SPECFEM3D_PMLs.txt
@@ -101,58 +101,6 @@ Thanks,
Dimitri.
------------------------------------------------
-
-Subject: Re: specfem3d (PS)
-Date: Sat, 24 Jan 2015 00:06:04 +0100
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: Zhang, Chang-Hua
-CC: Zhinan Xie
-
-Hi Chang-hua,
-
-PS: please note that more traditionally we tend to indicate "no
-attenuation" with a test such as Q > 9999. rather than Q < 0.0001 or Q == 0;
-this is because low values of Q mean high attenuation, thus no
-attenuation means Q tending to infinity rather than to zero.
-
-Thus, in the whole code tests such as
-
-if (Q_mu <=1.e-5) cycle
-
-should be changed to if(Q > 9999.) then
-
-and if the code detects Q = 0 (or Q < 1.e-5) in an input file it
-should instead stop with an error message, since a Q = 0 (or close to
-zero) value would be physically meaningless.
-
-Thanks,
-Best regards,
-
-Dimitri.
-
-------------------------------------------------
-
-Subject: Re: specfem3d
-Date: Sat, 24 Jan 2015 00:00:38 +0100
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: Zhang, Chang-Hua
-CC: Zhinan Xie
-
-Hi Chang-hua,
-
-Thank you very much! Your tests and your modifications are going to be
-very useful.
-
-Could you please send your modified routines to Zhinan, so that he can
-commit them to the official Git version of the code?
-
-Thank you,
-Best regards,
-
-Dimitri.
-
------------------------------------------------
Subject: RE: specfem3d
@@ -163,57 +111,14 @@ CC: Zhinan Xie
Hi Dimitri and Zhinan,
-Now I have done more tests, and it seems to me that by choosing Qs=35 in the example (homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides) is not appropriate because of the large element size, or the large resolved period. By creating a small element size (25 m) homogenous model, and Qs=35, the spike appearing in the snap shot shown in slide 5 disappears. So only the PML needs to be fix, which I have done that.
-
There are two more bugs:
1) The solver crashes for an acoustic-elastic coupling model with PML because of the reallocating an allocated arrays in pml_allocate_arrays_dummy(). In order to fix it, I have replaced the "call pml_allocate_arrays_dummy()" with if(.not.PML_CONDITIONS) call pml_allocate_arrays_dummy() in prepare_timerun.F90
-2) There is no Q_kappa input in the mesh_file for the internal mesher. Therefore, the "if (Q_kappa <=1.e-5) cycle" statement in get_attenuation_model.f90 always works, and not attenuation for both S and P waves even Q_mu is given. I think neither "if (Q_kappa <=1.e-5) cycle" nor "if (Q_mu <=1.e-5) cycle" serves the purpose if Q=0 means no attenuation.
-
-So my conclusion is that attenuation and CPML can now be both turned on.
Best Regards,
Chang-hua
------------------------------------------------
-
-Subject: Re: specfem3d
-Date: Fri, 23 Jan 2015 23:59:16 +0100
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: Zhang, Chang-Hua, xiezhinan
-
-Hi Chang-hua,
-
-Yes, you are right, in the current 3D code there is always this problem
-that attenuation and PML cannot be turned on simultaneously (Zhinan will
-soon remove that (unneeded, as you pointed out) constraint based on your
-suggestions.
-
-What Zhinan meant below in that in the 2D that unneeded constraint is
-not present and then when you use both viscoelasticity in the main
-domain and PML you get correct matching and correct absorption if Q is
-moderate (above 70 or so), however for very strong attenuation (Q below
-50 or so) then you start to see spurious reflections from the edge
-between the PML and the main domain because the main domain is
-viscoelastic but the PML is only elastic, thus there is a jump from Q =
-50 or so in the domain to Q = infinity in the PML, and that jump acts as
-an interface (an interface between two layers with different Q
-properties) that generates low-frequency reflections, as an absorbing
-sponge would do.
-
-Thus, in practice, using an elastic PML in contact with a viscoelastic
-domain leads to good results for moderate or low attenuation, but to
-results that start to show spurious reflections if attenuation is high
-
-(and that is normal; the only clean solution would be to implement PML
-for viscoelastic materials inside the PML, but as you said in your email
-that is complicated)
-
-Best regards,
-Dimitri.
-
-
------------------------------------------------
Subject: RE: Re:Re: specfem3d
@@ -236,186 +141,13 @@ Chang-hua
------------------------------------------------
-Subject: Re:Re: specfem3d
-Date: Fri, 23 Jan 2015 06:03:10 +0800 (CST)
-From: xiezhinan
-To: Dimitri Komatitsch
-CC: Zhang, Chang-Hua
-
-Dear Chang hua and Dimitri,
-
-Chang hua, first, thank you so much for your comment in CPML
- implementation in SPECFEM2D, which is not optimized.
-As you have already seen that we waste a lot of memory.
-If you do not mind, I will use your idea to reduce the memory cost in
-CPML implementation.
-
-I do not get your point of following paragraph
-
-Redefine the CPML memory variables: CPML memory variables for
-calculation of Txx, Txy, and Txzare labeled by _x, and for Tyx, Tyy, and
-Tyzlabeled by _y, and for Tzx, Tzy, and Tzzlabeled by _z. These
-redefinitions are convenient and necessary for the extension to the
-anisotropic case.
-
-Could you please make it more clear? I am new to anisotropic wave equations.
-
-Concerning the CPML in 3D do not work for viscoelastic wave simulation.
-We will fix that. But the way we current use is that we will use a PML
-formulated on basis of elastic wave equation to absorb the outgoing wave
-from interior domain governed by viscoelastic wave equation.
-
-Though the way is poor, it works when the Q is not too low. We have
-tested for Q>70, the result is quite well.
-
-Thank you so much.
-Best regards,
-Zhinan
-
------------------------------------------------
-
-Subject: Re: specfem3d
-Date: Thu, 22 Jan 2015 00:02:29 +0100
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: Zhang, Chang-Hua
-CC: Zhinan Xie
-
-Dear Chang-hua,
-
-Thanks! Let me forward this to Zhinan, so that he can have a look at
-page 4 of the PDF file you sent me, which shows that SPECFEM3D currently
-stops when one tries to use CPML when the interior of the domain has
-viscoelasticity.
-
-Thanks,
-Best regards,
-
-Dimitri.
-
--------------------------
-
-Subject: Re: [specfem3d] Runing both PML and attenuation causes error
-results (#356)
-Date: Tue, 20 Jan 2015 15:11:47 -0800
-From: specfem3d-zhang-ksu
-To: geodynamics/specfem3d
-CC: Dimitri Komatitsch
-
-Thanks, Dimitri. I have replaced the code segment by the following one,
-
-| ! if both PML and attenuation are used, the PML region must be elastic
- if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then
- ! do not merge this second line with the first using an ".and." statement
- ! because array is_CPML() is unallocated when PML_CONDITIONS is false
- if (is_CPML(ispec)) then
- ! PML region
- ispec_CPML = spec_to_CPML(ispec)
- DO_LOOP_IJK
-
dummyx_loc_att(INDEX_IJK) = PML_displ_old(1,INDEX_IJK,ispec_CPML)
dummyy_loc_att(INDEX_IJK) = PML_displ_old(2,INDEX_IJK,ispec_CPML)
dummyz_loc_att(INDEX_IJK) = PML_displ_old(3,INDEX_IJK,ispec_CPML)
- ENDDO_LOOP_IJK
- else if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
- ! attenuation region
- DO_LOOP_IJK
- iglob = ibool(INDEX_IJK,ispec)
- dummyx_loc_att(INDEX_IJK) = deltat*veloc(1,iglob)
- dummyy_loc_att(INDEX_IJK) = deltat*veloc(2,iglob)
- dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob)
-
- ENDDO_LOOP_IJK
- endif
- else if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
- ! only attenuation
- ! use first order Taylor expansion of displacement for local storage of stresses
- ! at this current time step, to fix attenuation in a consistent way
- DO_LOOP_IJK
- iglob = ibool(INDEX_IJK,ispec)
- dummyx_loc_att(INDEX_IJK) = deltat*veloc(1,iglob)
- dummyy_loc_att(INDEX_IJK) = deltat*veloc(2,iglob)
- dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob)
- ENDDO_LOOP_IJK
-
- endif
-|
-
-And similarly the other places, the results are much better, but I have
-a strong spike at the source location all the time, no idea what is
-going on.
-
----------------------------------------------------------
-Subject: Re: [specfem3d] Runing both PML and attenuation causes error results (#356)
-Date: Tue, 20 Jan 2015 23:58:21 +0100
-From: Dimitri Komatitsch
-Organization: CNRS, Marseille, France
-To: geodynamics/specfem3d , specfem3d-zhang-ksu
-CC: Zhinan Xie
-
-Hi,
-
-Yes, PML in the 3D code is currently for acoustic/elastic media only
-(i.e. viscoelasticity cannot be turned on).
-
-Let me cc my colleague Zhinan Xie because he is currently working on the
-code and thus he could probably remove that constraint (note that if so
-PML will *not* be perfectly matched strictly speaking because it is
-written without attenuation, thus the viscoelastic medium will be in
-contact with an elastic PML; but that works fine in practice if Q
-factors are not too small (typically not below 50 or so).
-
-PS for Zhinan: Zhang Ksu is right, there is no reason to stop the code
-in that case, we should just run with an elastic PML in contact with the
-viscoelastic medium, as we do in the 2D code. Thus the code in his email
-below should be slightly modified (and the PML arrays should be
-allocated even if attenuation is on I guess).
-
-Thanks,
-Dimitri.
-
-On 01/20/2015 12:00 AM, specfem3d-zhang-ksu wrote:
-> I am running the example of
-> homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides by turning
-> ATTENUATION false or true in the Par_file and change Qmu in the
-> nummaterial_velocity_file.
->
-> The result with ATTENUATION=.true. is totally wrong. Is it a problem
-> that the example is only setup for CPML or there is a bug in the code?
->
-> By examining the source code, I think that there is a bug is more
-> probable. For example, in the
-> subroutine compute_forces_viscoelastic_Dev_5p, there is code segment
-> ! use first order Taylor expansion of displacement for local storage of
-> stresses
-> ! at this current time step, to fix attenuation in a consistent way
-> if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
-> DO_LOOP_IJK
-> iglob = ibool(INDEX_IJK,ispec)
-> dummyx_loc_att(INDEX_IJK) = deltat/veloc(1,iglob)
-> dummyy_loc_att(INDEX_IJK) = deltat/veloc(2,iglob)
-> dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob)
-> ENDDO_LOOP_IJK
-> else if (PML_CONDITIONS .and. (.not. backward_simulation) .and.
-> NSPEC_CPML > 0) then
-> ! do not merge this second line with the first using an ".and." statement
-> ! because array is_CPML() is unallocated when PML_CONDITIONS is false
-> if (is_CPML(ispec)) then
-> DO_LOOP_IJK
-> iglob = ibool(INDEX_IJK,ispec)
-> dummyx_loc_att(INDEX_IJK) = displ_old(1,iglob)
-> dummyy_loc_att(INDEX_IJK) = displ_old(2,iglob)
-> dummyz_loc_att(INDEX_IJK) = displ_old(3,iglob)
-> ENDDO_LOOP_IJK
-> endif
-> endif
->
-> which seams to indicate we can not have both CPML and attenuation.
->
-> Has anyone experienced this issue?
-
==============================================================================
==============================================================================
==============================================================================
More information about the CIG-COMMITS
mailing list