[cig-commits] r18485 - in seismo/3D/CPML: tags/v1.1.3 trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat May 28 17:17:48 PDT 2011
Author: dkomati1
Date: 2011-05-28 17:17:47 -0700 (Sat, 28 May 2011)
New Revision: 18485
Modified:
seismo/3D/CPML/tags/v1.1.3/README_seismic_cpml.html
seismo/3D/CPML/tags/v1.1.3/seismic_ADEPML_2D_RK4_eighth_order.f90
seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_anisotropic.f90
seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_fourth_order.f90
seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_second_order.f90
seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_poroelastic_fourth_order.f90
seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_isotropic_MPI_OpenMP.f90
seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_viscoelastic_MPI.f90
seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_anisotropic_fourth.f90
seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_isotropic.f90
seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_3D_isotropic_OpenMP.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90
seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90
seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90
seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90
seismo/3D/CPML/trunk/seismic_PML_Collino_2D_anisotropic_fourth.f90
Log:
moved BibTeX reference KoMa07 further down in the list of references in the comments
Modified: seismo/3D/CPML/tags/v1.1.3/README_seismic_cpml.html
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/README_seismic_cpml.html 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/README_seismic_cpml.html 2011-05-29 00:17:47 UTC (rev 18485)
@@ -3,10 +3,10 @@
<HEAD>
<META HTTP-EQUIV="CONTENT-TYPE" CONTENT="text/html; charset=utf-8">
<TITLE>The SEISMIC_CPML software package</TITLE>
- <META NAME="GENERATOR" CONTENT="OpenOffice.org 3.1 (Linux)">
+ <META NAME="GENERATOR" CONTENT="OpenOffice.org 3.3 (Unix)">
<META NAME="CREATED" CONTENT="0;0">
<META NAME="CHANGEDBY" CONTENT="Dimitri Komatitsch">
- <META NAME="CHANGED" CONTENT="20101219;534200">
+ <META NAME="CHANGED" CONTENT="20110503;580000">
<STYLE TYPE="text/css">
<!--
P { color: #000000; font-family: "Times New Roman"; font-size: 12pt }
@@ -15,10 +15,10 @@
-->
</STYLE>
</HEAD>
-<BODY LANG="en-US" LINK="#0000ff" VLINK="#0000ff" BGCOLOR="#ffffff" BACKGROUND="http://www.univ-pau.fr/~dkomati1/grayback.gif" DIR="LTR">
-<P><A HREF="http://www.univ-pau.fr/~dkomati1">Home page of Dimitri
+<BODY LANG="en-US" LINK="#0000ff" VLINK="#0000ff" BGCOLOR="#ffffff" BACKGROUND="http://www.univ-pau.fr/%7Edkomati1/grayback.gif" DIR="LTR">
+<P><A HREF="http://www.univ-pau.fr/%7Edkomati1">Home page of Dimitri
Komatitsch</A></P>
-<P ALIGN=CENTER><A NAME="_x0000_i1025"></A><IMG SRC="http://www.univ-pau.fr/~dkomati1/seismic_cpml.gif" NAME="graphics1" ALIGN=BOTTOM WIDTH=158 HEIGHT=30 BORDER=0></P>
+<P ALIGN=CENTER><A NAME="_x0000_i1025"></A><IMG SRC="http://www.univ-pau.fr/%7Edkomati1/seismic_cpml.gif" NAME="graphics1" ALIGN=BOTTOM WIDTH=158 HEIGHT=30 BORDER=0></P>
<P> </P>
<P ALIGN=CENTER><A NAME="download_map"></A><IMG SRC="http://www.geodynamics.org/~buildbot/maps/Seismic_CPML.gif" NAME="download_map2" ALIGN=BOTTOM WIDTH=445 HEIGHT=223 BORDER=0></P>
<P ALIGN=CENTER><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><U><B>Current
@@ -28,13 +28,13 @@
is a set of ten open-source Fortran90 programs</FONT></FONT> <FONT SIZE=3 STYLE="font-size: 13pt">to
solve the two-dimensional or three-dimensional isotropic or
anisotropic elastic, viscoelastic or poroelastic wave equation using
-a finite-difference method with Convolutional Perfectly Matched Layer
-(C-PML) conditions, developed by Dimitri Komatitsch and Roland Martin
-from University of Pau, France. </FONT>
+a finite-difference method with Convolutional or Auxiliary Perfectly
+Matched Layer (C-PML or ADE-PML) conditions, developed by Dimitri
+Komatitsch and Roland Martin from University of Pau, France. </FONT>
</P>
<P STYLE="margin-bottom: 0in"><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt">You
can get the full source code of the programs at the official Web
-site: </FONT></FONT></FONT><A HREF="http://www.geodynamics.org/cig/software"><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><B>http://www.geodynamics.org/cig/software</B></FONT></FONT></FONT></A></P>
+site: </FONT></FONT></FONT><A HREF="http://www.geodynamics.org/cig/software"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt">http://www.geodynamics.org/cig/software</FONT></FONT></A></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">The unsplit <B>Convolutional
Perfectly Matched Layer (C-PML) for the 3D elastic wave equation</B>
was introduced and is described in detail in: </FONT>
@@ -44,8 +44,8 @@
unsplit convolutional Perfectly Matched Layer improved at grazing
incidence for the seismic wave equation</SPAN></FONT><FONT SIZE=3 STYLE="font-size: 13pt"><I>,
Geophysics</I></FONT><FONT SIZE=3 STYLE="font-size: 13pt">, vol.
-72(5), p SM155-SM167, doi: 10.1190/1.2757586 (2007). <A HREF="http://www.univ-pau.fr/~dkomati1/published_papers/geophysics_CPML_2007_elastic_typos_fixed.pdf">PDF
-reprint</A> <A HREF="http://www.univ-pau.fr/~dkomati1/bibtex_komatitsch.bib">BibTeX</A></FONT></P>
+72(5), p SM155-SM167, doi: 10.1190/1.2757586 (2007). <A HREF="http://www.univ-pau.fr/%7Edkomati1/published_papers/geophysics_CPML_2007_elastic_typos_fixed.pdf">PDF
+reprint</A> <A HREF="http://www.univ-pau.fr/%7Edkomati1/bibtex_komatitsch.bib">BibTeX</A></FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">It was originally developed
for Maxwell's equations by Roden and Gedney (2000) (see reference
below).</FONT></P>
@@ -59,8 +59,8 @@
</SPAN></FONT></SPAN><SPAN STYLE="text-decoration: none"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN LANG="fr-FR"><SPAN STYLE="font-weight: normal">vol.
179(1), p. 333-344, </SPAN></SPAN></FONT></SPAN><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN LANG="fr-FR"><SPAN STYLE="font-style: normal"><SPAN STYLE="font-weight: normal">doi:
10.1111/j.1365-246X.2009.04278.x </SPAN></SPAN></SPAN></FONT></FONT></SPAN></FONT><SPAN STYLE="text-decoration: none"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN STYLE="font-weight: normal">(2009).</SPAN></FONT></SPAN><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none">
-</SPAN></FONT><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN LANG="fr-FR"><SPAN STYLE="font-weight: normal"><A HREF="http://www.univ-pau.fr/~dkomati1/published_papers/GJI_CPML_2009_viscoelastic.pdf">PDF
-reprint</A> <A HREF="http://www.univ-pau.fr/~dkomati1/bibtex_komatitsch.bib">BibTeX</A></SPAN></SPAN></FONT></FONT></SPAN></FONT></P>
+</SPAN></FONT><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN LANG="fr-FR"><SPAN STYLE="font-weight: normal"><A HREF="http://www.univ-pau.fr/%7Edkomati1/published_papers/GJI_CPML_2009_viscoelastic.pdf">PDF
+reprint</A> <A HREF="http://www.univ-pau.fr/%7Edkomati1/bibtex_komatitsch.bib">BibTeX</A></SPAN></SPAN></FONT></FONT></SPAN></FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">An extension to poroelastic
media is developed in:</FONT></P>
<P><SPAN STYLE="text-decoration: none"><FONT SIZE=3 STYLE="font-size: 13pt"><B>Roland
@@ -69,8 +69,8 @@
unsplit convolutional Perfectly Matched Layer improved at grazing
incidence for seismic wave propagation in poroelastic media</SPAN></FONT><FONT SIZE=3 STYLE="font-size: 13pt"><I>,
Geophysics</I></FONT><FONT SIZE=3 STYLE="font-size: 13pt">, vol.
-73(4), p T51-T61, doi: 10.1190/1.2939484 (2008). <A HREF="http://www.univ-pau.fr/~dkomati1/published_papers/geophysics_CPML_2008_poroelastic_typos_fixed.pdf">PDF
-reprint</A> <A HREF="http://www.univ-pau.fr/~dkomati1/bibtex_komatitsch.bib">BibTeX</A></FONT></P>
+73(4), p T51-T61, doi: 10.1190/1.2939484 (2008). <A HREF="http://www.univ-pau.fr/%7Edkomati1/published_papers/geophysics_CPML_2008_poroelastic_typos_fixed.pdf">PDF
+reprint</A> <A HREF="http://www.univ-pau.fr/%7Edkomati1/bibtex_komatitsch.bib">BibTeX</A></FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">and a variational formulation
is developed in:</FONT></P>
<P><SPAN STYLE="text-decoration: none"><FONT SIZE=3 STYLE="font-size: 13pt"><B>Roland
@@ -82,8 +82,8 @@
perfectly matched layer for the isotropic or anisotropic seismic wave
equation, </FONT><FONT SIZE=3 STYLE="font-size: 13pt"><I>Computer
Modeling in Engineering and Sciences</I></FONT><FONT SIZE=3 STYLE="font-size: 13pt">,
-vol. 37(3), p. 274-304 (2008). </FONT><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN STYLE="font-weight: normal"><A HREF="http://www.univ-pau.fr/~dkomati1/published_papers/CMES_cpml_2008.pdf">PDF
-reprint</A> <A HREF="http://www.univ-pau.fr/~dkomati1/bibtex_komatitsch.bib">BibTeX</A></SPAN></FONT></FONT></SPAN></FONT></P>
+vol. 37(3), p. 274-304 (2008). </FONT><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=3 STYLE="font-size: 13pt"><SPAN STYLE="font-weight: normal"><A HREF="http://www.univ-pau.fr/%7Edkomati1/published_papers/CMES_cpml_2008.pdf">PDF
+reprint</A> <A HREF="http://www.univ-pau.fr/%7Edkomati1/bibtex_komatitsch.bib">BibTeX</A></SPAN></FONT></FONT></SPAN></FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">An extension to higher-order
time schemes, called ADE-PML (Auxiliary Differential Equation - PML)
is developed in:</FONT></P>
@@ -93,9 +93,9 @@
matched layer for the seismic wave equation using Auxiliary
Differential Equations (ADE-PML), </SPAN></SPAN></FONT></FONT></SPAN></FONT><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=4><SPAN LANG="fr-FR"><I><SPAN STYLE="font-weight: normal">Computer
Modeling in Engineering and Sciences</SPAN></I></SPAN></FONT></FONT></SPAN></FONT><FONT COLOR="#000000"><SPAN STYLE="text-decoration: none"><FONT FACE="Times New Roman, serif"><FONT SIZE=4><SPAN LANG="fr-FR"><SPAN STYLE="font-weight: normal">,
-vol. 56(1), p. 17-42 (2010).</SPAN></SPAN></FONT></FONT></SPAN></FONT><FONT SIZE=3 STYLE="font-size: 13pt">
-<A HREF="http://www.univ-pau.fr/~dkomati1/published_papers/CMES_ADE_PML_2010.pdf">PDF
-reprint</A> <A HREF="http://www.univ-pau.fr/~dkomati1/bibtex_komatitsch.bib">BibTeX</A></FONT></P>
+vol. 56(1), p. 17-42 (2010).</SPAN></SPAN></FONT></FONT></SPAN></FONT>
+<FONT SIZE=3 STYLE="font-size: 13pt"><A HREF="http://www.univ-pau.fr/%7Edkomati1/published_papers/CMES_ADE_PML_2010.pdf">PDF
+reprint</A> <A HREF="http://www.univ-pau.fr/%7Edkomati1/bibtex_komatitsch.bib">BibTeX</A></FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">Note that in the case of an
anisotropic medium the modification made is not strictly speaking
perfectly matched any more, i.e., not a PML, but rather </FONT><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt">a
@@ -105,12 +105,12 @@
vol. 98(4), p. 1811-1836 (2008). H</FONT></FONT></FONT><FONT SIZE=3 STYLE="font-size: 13pt">owever,
it works well in practice even if it is not perfectly matched any
more from a mathematical point of view.</FONT></P>
-<P><FONT COLOR="#ff0000"><B><FONT SIZE=3 STYLE="font-size: 13pt">IMPORTANT:
+<P><FONT COLOR="#ff0000"><FONT SIZE=3 STYLE="font-size: 13pt"><B>IMPORTANT:
all of our codes are written in Fortran; if you have written or if
you write a C or C++ version of some of these codes and want to make
them open source (GNU GPL v2 or CeCILL, as you prefer) and part of
the package, please do not hesitate to send them to us, we will add
-them to our tar file and will acknowledge you as the author.</FONT></B></FONT></P>
+them to our tar file and will acknowledge you as the author.</B></FONT></FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">This software is governed by
the <A HREF="http://www.cecill.info/licences/Licence_CeCILL_V2-en.html">CeCILL
license (version 2)</A> (a French license very similar to GNU GPL
@@ -120,18 +120,12 @@
CNRS and INRIA at the following URL "<A HREF="http://www.cecill.info/index.en.html">http://www.cecill.info</A>".</FONT></P>
<P><FONT SIZE=3 STYLE="font-size: 13pt">If you use this code for your
own research, please cite some (or all) of these articles:</FONT></P>
-<P STYLE="margin-bottom: 0in"><FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">@ARTICLE{KoMa07,<BR>author
-= {Dimitri Komatitsch and Roland Martin},<BR>title = {An unsplit
-convolutional {P}erfectly {M}atched {L}ayer improved at grazing
-incidence for the seismic wave equation},<BR>journal =
-{Geophysics},<BR>year = {2007},<BR>volume = {72},<BR>number =
-{5},<BR>pages = {SM155-SM167},<BR>doi =
-{10.1190/1.2757586}}<BR><BR>@ARTICLE{MaKo09,<BR>author = {Roland
-Martin and Dimitri Komatitsch},<BR>title = {An unsplit convolutional
-perfectly matched layer technique improved at grazing incidence for
-the viscoelastic wave equation},<BR>journal = {Geophysical Journal
-International},<BR>year = {2009},<BR>volume = {179},<BR>number =
-{1},<BR>pages = {333-344},<BR>doi =
+<P><FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">@ARTICLE{MaKo09,<BR>author
+= {Roland Martin and Dimitri Komatitsch},<BR>title = {An unsplit
+convolutional perfectly matched layer technique improved at grazing
+incidence for the viscoelastic wave equation},<BR>journal =
+{Geophysical Journal International},<BR>year = {2009},<BR>volume =
+{179},<BR>number = {1},<BR>pages = {333-344},<BR>doi =
{10.1111/j.1365-246X.2009.04278.x}}<BR><BR>@ARTICLE{MaKoEz08,<BR>author
= {Roland Martin and Dimitri Komatitsch and Abdelaaziz
Ezziani},<BR>title = {An unsplit convolutional perfectly matched
@@ -145,37 +139,21 @@
equation},<BR>journal = {Computer Modeling in Engineering and
Sciences},<BR>year = {2008},<BR>volume = {37},<BR>pages =
{274-304},<BR>number = {3}}</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"><BR>
-</P>
-<P STYLE="margin-bottom: 0in"><FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">@ARTICLE{MaKoGeBr10,
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">author
+<P STYLE="margin-bottom: 0in"><FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">@ARTICLE{MaKoGeBr10,<BR>author
= {Roland Martin and Dimitri Komatitsch and Stephen D. Gedney and
-Emilien
- Bruthiaux},
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">title
-= {A high-order time and space formulation of the unsplit perfectly
-matched layer for the seismic wave equation using {Auxiliary
-Differential
- Equations (ADE-PML)}},
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">journal
-= {Computer Modeling in Engineering and Sciences},
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">year
-= {2010},
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">volume
-= {56},
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">pages
-= {17-42},
-</FONT></FONT></P>
-<P STYLE="margin-bottom: 0in"> <FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">number
-= {1}
-}
-<BR><BR></FONT></FONT><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt">Roden
+Emilien Bruthiaux},<BR>title = {A high-order time and space
+formulation of the unsplit perfectly matched layer for the seismic
+wave equation using {Auxiliary Differential Equations
+(ADE-PML)}},<BR>journal = {Computer Modeling in Engineering and
+Sciences},<BR>year = {2010},<BR>volume = {56},<BR>pages =
+{17-42},<BR>number = {1}}</FONT></FONT></P>
+<P STYLE="margin-bottom: 0in"><FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">@ARTICLE{KoMa07,<BR>author
+= {Dimitri Komatitsch and Roland Martin},<BR>title = {An unsplit
+convolutional {P}erfectly {M}atched {L}ayer improved at grazing
+incidence for the seismic wave equation},<BR>journal =
+{Geophysics},<BR>year = {2007},<BR>volume = {72},<BR>number =
+{5},<BR>pages = {SM155-SM167},<BR>doi = {10.1190/1.2757586}}</FONT></FONT></P>
+<P STYLE="margin-bottom: 0in"><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=3 STYLE="font-size: 13pt">Roden
and Gedney's original article for Maxwell's equations
is:</FONT></FONT></FONT><FONT COLOR="#000000"><FONT FACE="Liberation Serif, serif"><FONT SIZE=1 STYLE="font-size: 6pt"><BR><BR><BR></FONT></FONT></FONT><FONT FACE="Courier New, monospace"><FONT SIZE=3 STYLE="font-size: 13pt">@ARTICLE{RoGe00,<BR>author
= {J. A. Roden and S. D. Gedney},<BR>title = {Convolution {PML}
@@ -241,7 +219,7 @@
text of the CeCILL licence (version 2)</FONT></A></P>
<P><A HREF="http://www.cecill.info/index.en.html"><FONT SIZE=3 STYLE="font-size: 13pt">Official
web site of the CeCILL licence</FONT></A></P>
-<P><A HREF="http://www.univ-pau.fr/~dkomati1"><FONT SIZE=3 STYLE="font-size: 13pt">Home
+<P><A HREF="http://www.univ-pau.fr/%7Edkomati1"><FONT SIZE=3 STYLE="font-size: 13pt">Home
page of Dimitri Komatitsch</FONT></A></P>
</BODY>
-</HTML>
\ No newline at end of file
+</HTML>
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_ADEPML_2D_RK4_eighth_order.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_ADEPML_2D_RK4_eighth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_ADEPML_2D_RK4_eighth_order.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -327,16 +327,16 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(-4:NX+4) :: d_x_1,K_x_1,alpha_prime_x_1,g_x_1,ksi_x
- double precision, dimension(-4:NX+4) :: d_x_half_1,K_x_half_1,alpha_prime_x_half_1,g_x_half_1,ksi_x_half
- double precision, dimension(-4:NY+4) :: d_y_1,K_y_1,alpha_prime_y_1,g_y_1,ksi_y
- double precision, dimension(-4:NY+4) :: d_y_half_1,K_y_half_1,alpha_prime_y_half_1,g_y_half_1,ksi_y_half
+ double precision, dimension(-4:NX+4) :: d_x_1,K_x_1,alpha_x_1,g_x_1,ksi_x
+ double precision, dimension(-4:NX+4) :: d_x_half_1,K_x_half_1,alpha_x_half_1,g_x_half_1,ksi_x_half
+ double precision, dimension(-4:NY+4) :: d_y_1,K_y_1,alpha_y_1,g_y_1,ksi_y
+ double precision, dimension(-4:NY+4) :: d_y_half_1,K_y_half_1,alpha_y_half_1,g_y_half_1,ksi_y_half
! 1D arrays for the damping profiles
- double precision, dimension(-4:NX+4) :: d_x_2,K_x_2,alpha_prime_x_2,g_x_2
- double precision, dimension(-4:NX+4) :: d_x_half_2,K_x_half_2,alpha_prime_x_half_2,g_x_half_2
- double precision, dimension(-4:NY+4) :: d_y_2,K_y_2,alpha_prime_y_2,g_y_2
- double precision, dimension(-4:NY+4) :: d_y_half_2,K_y_half_2,alpha_prime_y_half_2,g_y_half_2
+ double precision, dimension(-4:NX+4) :: d_x_2,K_x_2,alpha_x_2,g_x_2
+ double precision, dimension(-4:NX+4) :: d_x_half_2,K_x_half_2,alpha_x_half_2,g_x_half_2
+ double precision, dimension(-4:NY+4) :: d_y_2,K_y_2,alpha_y_2,g_y_2
+ double precision, dimension(-4:NY+4) :: d_y_half_2,K_y_half_2,alpha_y_half_2,g_y_half_2
! coefficients that allow to reset the memory variables at each RK4 substep depend on the substepping and are then of dimension 4,
! 1D arrays for the damping profiles
@@ -409,13 +409,13 @@
thickness_PML_x = NPOINTS_PML * DELTAX
thickness_PML_y = NPOINTS_PML * DELTAY
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.00001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0_x = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_y)
@@ -440,8 +440,8 @@
d_x_half_1(:) = ZERO
K_x_1(:) = 1.d0
K_x_half_1(:) = 1.d0
- alpha_prime_x_1(:) = ZERO
- alpha_prime_x_half_1(:) = ZERO
+ alpha_x_1(:) = ZERO
+ alpha_x_half_1(:) = ZERO
a_x_1(:,:) = ZERO
a_x_half_1(:,:) = ZERO
g_x_1(:) = 5.d-1
@@ -453,8 +453,8 @@
d_y_half_1(:) = ZERO
K_y_1(:) = 1.d0
K_y_half_1(:) = 1.d0
- alpha_prime_y_1(:) = ZERO
- alpha_prime_y_half_1(:) = ZERO
+ alpha_y_1(:) = ZERO
+ alpha_y_half_1(:) = ZERO
a_y_1(:,:) = ZERO
a_y_half_1(:,:) = ZERO
g_y_1(:) = 1.d0
@@ -464,8 +464,8 @@
d_x_half_2(:) = ZERO
K_x_2(:) = 1.d0
K_x_half_2(:) = 1.d0
- alpha_prime_x_2(:) = ZERO
- alpha_prime_x_half_2(:) = ZERO
+ alpha_x_2(:) = ZERO
+ alpha_x_half_2(:) = ZERO
a_x_2(:,:) = ZERO
a_x_half_2(:,:) = ZERO
g_x_2(:) = 1.d0
@@ -475,8 +475,8 @@
d_y_half_2(:) = ZERO
K_y_2(:) = 1.d0
K_y_half_2(:) = 1.d0
- alpha_prime_y_2(:) = ZERO
- alpha_prime_y_half_2(:) = ZERO
+ alpha_y_2(:) = ZERO
+ alpha_y_half_2(:) = ZERO
a_y_2(:,:) = ZERO
a_y_half_2(:,:) = ZERO
g_y_2(:) = 1.d0
@@ -503,7 +503,7 @@
d_x_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -513,7 +513,7 @@
d_x_half_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
@@ -528,7 +528,7 @@
d_x_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -538,7 +538,7 @@
d_x_half_1(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_1(i) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_x_half_1(i) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
@@ -548,29 +548,29 @@
d_x_half_2(i) = 0.d0
! just in case, for -5 at the end
- if(alpha_prime_x_1(i) < ZERO) alpha_prime_x_1(i) = ZERO
- if(alpha_prime_x_half_1(i) < ZERO) alpha_prime_x_half_1(i) = ZERO
+ if(alpha_x_1(i) < ZERO) alpha_x_1(i) = ZERO
+ if(alpha_x_half_1(i) < ZERO) alpha_x_half_1(i) = ZERO
! just in case, for -5 at the end
- if(alpha_prime_x_2(i) < ZERO) alpha_prime_x_2(i) = ZERO
- if(alpha_prime_x_half_2(i) < ZERO) alpha_prime_x_half_2(i) = ZERO
+ if(alpha_x_2(i) < ZERO) alpha_x_2(i) = ZERO
+ if(alpha_x_half_2(i) < ZERO) alpha_x_half_2(i) = ZERO
! CPML damping parameters for the 4 sub time steps of RK4 algorithm
do inc=1,4
- b_x_1(inc,i) = (1.-epsn*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_prime_x_1(i)))/&
- (1.+epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_prime_x_1(i)))
+ b_x_1(inc,i) = (1.-epsn*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_x_1(i)))/&
+ (1.+epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_x_1(i)))
b_x_half_1(inc,i) = (1.-epsn*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i) &
- + alpha_prime_x_half_1(i)))/(1. +epsn1*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i) &
- + alpha_prime_x_half_1(i)))
+ + alpha_x_half_1(i)))/(1. +epsn1*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i) &
+ + alpha_x_half_1(i)))
! this to avoid division by zero outside the PML
if(abs(d_x_1(i)) > 1.d-6) a_x_1(inc,i) = - DELTAT*rk41(inc)*d_x_1(i) / (K_x_1(i)* K_x_1(i))/&
- (1. +epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_prime_x_1(i)))
+ (1. +epsn1*DELTAT*rk41(inc)*(d_x_1(i)/K_x_1(i) + alpha_x_1(i)))
if(abs(d_x_half_1(i)) > 1.d-6) a_x_half_1(inc,i) =-DELTAT*rk41(inc)*d_x_half_1(i)/&
(K_x_half_1(i)*K_x_half_1(i) )/&
(1. +epsn1*DELTAT*rk41(inc)*(d_x_half_1(i)/K_x_half_1(i)&
- + alpha_prime_x_half_1(i)))
+ + alpha_x_half_1(i)))
enddo
@@ -597,7 +597,7 @@
d_y_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -607,7 +607,7 @@
d_y_half_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
@@ -622,7 +622,7 @@
d_y_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -632,30 +632,30 @@
d_y_half_1(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_1(j) = 1.d0 + (K_MAX_PML_1 - 1.d0) * abscissa_normalized**NPOWER2
- alpha_prime_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
+ alpha_y_half_1(j) = ALPHA_MAX_PML_1 * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_y_1(j) < ZERO) alpha_prime_y_1(j) = ZERO
- if(alpha_prime_y_half_1(j) < ZERO) alpha_prime_y_half_1(j) = ZERO
+ if(alpha_y_1(j) < ZERO) alpha_y_1(j) = ZERO
+ if(alpha_y_half_1(j) < ZERO) alpha_y_half_1(j) = ZERO
! CPML damping parameters for the 4 sub time steps of RK4 algorithm
do inc=1,4
- b_y_1(inc,j) = (1.-epsn*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_prime_y_1(j)))/&
- (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_prime_y_1(j)))
+ b_y_1(inc,j) = (1.-epsn*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_y_1(j)))/&
+ (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_y_1(j)))
b_y_half_1(inc,j) = (1.-epsn*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j) + &
- alpha_prime_y_half_1(j)))/(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j)&
- + alpha_prime_y_half_1(j)))
+ alpha_y_half_1(j)))/(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j)&
+ + alpha_y_half_1(j)))
! this to avoid division by zero outside the PML
if(abs(d_y_1(j)) > 1.d-6) a_y_1(inc,j) = - DELTAT*rk41(inc)*d_y_1(j)&
/ (K_y_1(j)* K_y_1(j))/&
- (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_prime_y_1(j)))
+ (1.+epsn1*DELTAT*rk41(inc)*(d_y_1(j)/K_y_1(j) + alpha_y_1(j)))
if(abs(d_y_half_1(j)) > 1.d-6) a_y_half_1(inc,j) = -DELTAT*rk41(inc)*d_y_half_1(j) /&
(K_y_half_1(j) * K_y_half_1(j) )/&
-(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j) + alpha_prime_y_half_1(j)))
+(1.+epsn1*DELTAT*rk41(inc)*(d_y_half_1(j)/K_y_half_1(j) + alpha_y_half_1(j)))
enddo
enddo
@@ -1045,9 +1045,9 @@
! check stability of the code, exit if unstable
if(velocnorm > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
- call create_2D_image(vx(1:NX,1:NY),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx(1:NX,1:NY),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,1)
- call create_2D_image(vy(1:NX,1:NY),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy(1:NX,1:NY),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,2)
open(unit=20,file='energy.dat',status='unknown')
do it2 = 1,NSTEP
@@ -1185,7 +1185,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -1218,6 +1218,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -1322,7 +1323,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_anisotropic.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_anisotropic.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_anisotropic.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -71,17 +71,6 @@
! If you use this code for your own research, please cite some (or all) of these
! articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoGe08,
! author = {Roland Martin and Dimitri Komatitsch and Stephen D. Gedney},
! title = {A variational formulation of a stabilized unsplit convolutional perfectly
@@ -114,6 +103,17 @@
! number = {1},
! doi = {10.1111/j.1365-246X.2009.04278.x}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! If you use the anisotropic implementation, please cite this article,
! in which the anisotropic parameters are described, as well:
!
@@ -279,7 +279,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -304,8 +304,8 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
double precision :: thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision :: Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -346,7 +346,7 @@
print *,'Velocity of qSV along horizontal axis . . =',sqrt(c33/rho)
print *
-! from Becache et al., INRIA report, equation 7 page 5
+! from Becache et al., INRIA report, equation 7 page 5 http://hal.inria.fr/docs/00/07/22/83/PDF/RR-4304.pdf
if(c11*c22 - c12*c12 <= 0.d0) stop 'problem in definition of orthotropic material'
! check intrinsic mathematical stability of PML model for an anisotropic material
@@ -379,13 +379,13 @@
thickness_PML_x = NPOINTS_PML * DELTAX
thickness_PML_y = NPOINTS_PML * DELTAY
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0_x = - (NPOWER + 1) * quasi_cp_max * log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * quasi_cp_max * log(Rcoef) / (2.d0 * thickness_PML_y)
@@ -397,8 +397,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -406,8 +406,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -432,7 +432,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -442,7 +442,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -457,7 +457,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -467,22 +467,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -507,7 +507,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -517,7 +517,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -532,7 +532,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -542,18 +542,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -722,9 +722,9 @@
! check stability of the code, exit if unstable
if(velocnorm > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
- call create_2D_image(vx,NX,NY,it,ISOURCE,JSOURCE, &
+ call create_color_image(vx,NX,NY,it,ISOURCE,JSOURCE, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,1)
- call create_2D_image(vy,NX,NY,it,ISOURCE,JSOURCE, &
+ call create_color_image(vy,NX,NY,it,ISOURCE,JSOURCE, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,2)
endif
@@ -742,7 +742,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -773,6 +773,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -864,7 +865,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_fourth_order.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_fourth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_fourth_order.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -77,17 +77,6 @@
!
! If you use this code for your own research, please cite some (or all) of these articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdelaaziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -120,6 +109,17 @@
! pages = {334-339},
! doi = {10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! To display the 2D results as color images, use:
!
! " display image*.gif " or " gimp image*.gif "
@@ -221,7 +221,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -246,8 +246,8 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
double precision :: thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision :: Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -292,13 +292,13 @@
thickness_PML_x = NPOINTS_PML * DELTAX
thickness_PML_y = NPOINTS_PML * DELTAY
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0_x = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_y)
@@ -310,8 +310,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -319,8 +319,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -345,7 +345,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -355,7 +355,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -370,7 +370,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -380,22 +380,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -420,7 +420,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -430,7 +430,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -445,7 +445,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -455,18 +455,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -731,9 +731,9 @@
! check stability of the code, exit if unstable
if(velocnorm > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
- call create_2D_image(vx,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,1)
- call create_2D_image(vy,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,2)
endif
@@ -864,7 +864,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -897,6 +897,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -1001,7 +1002,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_second_order.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_second_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_isotropic_second_order.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -69,17 +69,6 @@
! If you use this code for your own research, please cite some (or all) of these
! articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdela\^aziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -112,6 +101,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
@@ -220,7 +220,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -245,8 +245,8 @@
value_dsigmaxy_dy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
double precision :: thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision :: Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -291,13 +291,13 @@
thickness_PML_x = NPOINTS_PML * DELTAX
thickness_PML_y = NPOINTS_PML * DELTAY
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0_x = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_y)
@@ -309,8 +309,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -318,8 +318,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -344,7 +344,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -354,7 +354,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -369,7 +369,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -379,22 +379,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -419,7 +419,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -429,7 +429,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -444,7 +444,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -454,18 +454,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -730,9 +730,9 @@
! check stability of the code, exit if unstable
if(velocnorm > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
- call create_2D_image(vx,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,1)
- call create_2D_image(vy,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,2)
endif
@@ -863,7 +863,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -896,6 +896,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -1000,7 +1001,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_poroelastic_fourth_order.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_poroelastic_fourth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_CPML_2D_poroelastic_fourth_order.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -82,17 +82,6 @@
! number = {4},
! doi = {10.1190/1.2939484}}
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKo09,
! author = {Roland Martin and Dimitri Komatitsch},
! title = {An unsplit convolutional perfectly matched layer technique improved
@@ -114,6 +103,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
@@ -292,7 +292,7 @@
! double precision, parameter :: K_MAX_PML = 7.d0 ! from Gedney page 8.11
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
! double precision, parameter :: ALPHA_MAX_PML = 0.05d0 ! from Gedney page 8.22
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! 2D arrays for the memory variables
double precision, dimension(0:NX+1,0:NY+1) :: gamma11,gamma22
@@ -306,8 +306,8 @@
memory_dy_sigmaxx,memory_dy_sigmayy,memory_dy_sigmaxy
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half_x,K_x_half_x,alpha_prime_x_half_x,a_x_half_x,b_x_half_x
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half_y,K_y_half_y,alpha_prime_y_half_y,a_y_half_y,b_y_half_y
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half_x,K_x_half_x,alpha_x_half_x,a_x_half_x,b_x_half_x
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half_y,K_y_half_y,alpha_y_half_y,a_y_half_y,b_y_half_y
double precision thickness_PML_x,thickness_PML_y,xoriginleft,xoriginright,yoriginbottom,yorigintop
double precision Rcoef,d0_x,d0_y,xval,yval,abscissa_in_PML,abscissa_normalized
@@ -391,13 +391,13 @@
thickness_PML_x = NPOINTS_PML * DELTAX
thickness_PML_y = NPOINTS_PML * DELTAY
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
if(HETEROGENEOUS_MODEL) then
d0_x = - (NPOWER + 1) * max(cp_bottom,cp_top) * log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * max(cp_bottom,cp_top) * log(Rcoef) / (2.d0 * thickness_PML_y)
@@ -421,11 +421,11 @@
K_y(:) = 1.d0
K_y_half_y(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half_x(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half_x(:) = ZERO
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half_y(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half_y(:) = ZERO
a_x(:) = ZERO
a_x_half_x(:) = ZERO
@@ -452,7 +452,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -462,7 +462,7 @@
d_x_half_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -477,7 +477,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -487,23 +487,23 @@
d_x_half_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half_x(i) < ZERO) alpha_prime_x_half_x(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half_x(i) < ZERO) alpha_x_half_x(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half_x(i) = exp(- (d_x_half_x(i) / K_x_half_x(i) + alpha_prime_x_half_x(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half_x(i) = exp(- (d_x_half_x(i) / K_x_half_x(i) + alpha_x_half_x(i)) * DELTAT)
! this to avoid division by zero outside the PML
if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) /&
- (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half_x(i)) > 1.d-6) a_x_half_x(i) = d_x_half_x(i)&
- * (b_x_half_x(i) - 1.d0) / (K_x_half_x(i) * (d_x_half_x(i) + K_x_half_x(i) * alpha_prime_x_half_x(i)))
+ * (b_x_half_x(i) - 1.d0) / (K_x_half_x(i) * (d_x_half_x(i) + K_x_half_x(i) * alpha_x_half_x(i)))
enddo
@@ -528,7 +528,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -538,7 +538,7 @@
d_y_half_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -553,7 +553,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -563,23 +563,23 @@
d_y_half_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
-! if(alpha_prime_y(j) < ZERO) alpha_prime_y(j) = ZERO
-! if(alpha_prime_y_half_y(j) < ZERO) alpha_prime_y_half_y(j) = ZERO
+! if(alpha_y(j) < ZERO) alpha_y(j) = ZERO
+! if(alpha_y_half_y(j) < ZERO) alpha_y_half_y(j) = ZERO
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half_y(j) = exp(- (d_y_half_y(j) / K_y_half_y(j) + alpha_prime_y_half_y(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half_y(j) = exp(- (d_y_half_y(j) / K_y_half_y(j) + alpha_y_half_y(j)) * DELTAT)
! this to avoid division by zero outside the PML
if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) &
- / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half_y(j)) > 1.d-6) a_y_half_y(j) = d_y_half_y(j)&
- * (b_y_half_y(j) - 1.d0) / (K_y_half_y(j) * (d_y_half_y(j) + K_y_half_y(j) * alpha_prime_y_half_y(j)))
+ * (b_y_half_y(j) - 1.d0) / (K_y_half_y(j) * (d_y_half_y(j) + K_y_half_y(j) * alpha_y_half_y(j)))
enddo
@@ -969,11 +969,11 @@
vnorm(:,:)=sqrt(vx(:,:)**2+vy(:,:)**2)
- call create_2D_image(vx,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_LEFT,USE_PML_RIGHT,USE_PML_BOTTOM,&
USE_PML_TOP,1,max_amplitude,JINTERFACE)
- call create_2D_image(vy,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy,NX+2,NY+2,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_LEFT,USE_PML_RIGHT,USE_PML_BOTTOM,&
USE_PML_TOP,2,max_amplitude,JINTERFACE)
@@ -1111,7 +1111,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_LEFT,USE_PML_RIGHT,USE_PML_BOTTOM,USE_PML_TOP,field_number,max_amplitude,JINTERFACE)
@@ -1146,6 +1146,7 @@
integer :: R, G, B
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i5.5,'_Vx.pnm')") it
write(system_command,"('convert image',i5.5,'_Vx.pnm image',i5.5,'_Vx.gif ; rm image',i5.5,'_Vx.pnm')") it,it,it
@@ -1258,7 +1259,7 @@
! call the system to convert image to JPEG
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_isotropic_MPI_OpenMP.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_isotropic_MPI_OpenMP.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_isotropic_MPI_OpenMP.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -53,17 +53,6 @@
! If you use this code for your own research, please cite some (or all) of these
! articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdela\^aziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -96,6 +85,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
@@ -208,7 +208,7 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 1.d0 ! from Gedney page 8.11
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -253,9 +253,9 @@
value_dsigmayz_dz
! 1D arrays for the damping profiles
- double precision, dimension(NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
- double precision, dimension(NZ) :: d_z,K_z,alpha_prime_z,a_z,b_z,d_z_half,K_z_half,alpha_prime_z_half,a_z_half,b_z_half
+ double precision, dimension(NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
+ double precision, dimension(NZ) :: d_z,K_z,alpha_z,a_z,b_z,d_z_half,K_z_half,alpha_z_half,a_z_half,b_z_half
! PML
double precision thickness_PML_x,thickness_PML_y,thickness_PML_z
@@ -396,13 +396,13 @@
thickness_PML_y = NPOINTS_PML * DELTAY
thickness_PML_z = NPOINTS_PML * DELTAZ
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0_x = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_y)
d0_z = - (NPOWER + 1) * cp * log(Rcoef) / (2.d0 * thickness_PML_z)
@@ -419,8 +419,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -428,8 +428,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -437,8 +437,8 @@
d_z_half(:) = ZERO
K_z(:) = 1.d0
K_z_half(:) = 1.d0
- alpha_prime_z(:) = ZERO
- alpha_prime_z_half(:) = ZERO
+ alpha_z(:) = ZERO
+ alpha_z_half(:) = ZERO
a_z(:) = ZERO
a_z_half(:) = ZERO
@@ -463,7 +463,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -473,7 +473,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -488,7 +488,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -498,22 +498,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -538,7 +538,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -548,7 +548,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -563,7 +563,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -573,18 +573,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -609,7 +609,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -619,7 +619,7 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -634,7 +634,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -644,18 +644,18 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_prime_z(k)) * DELTAT)
- b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_prime_z_half(k)) * DELTAT)
+ b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_z(k)) * DELTAT)
+ b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_z_half(k)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_prime_z(k)))
+ if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_z(k)))
if(abs(d_z_half(k)) > 1.d-6) a_z_half(k) = d_z_half(k) * &
- (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_prime_z_half(k)))
+ (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_z_half(k)))
enddo
@@ -1223,9 +1223,9 @@
print *
call write_seismograms(sisvx,sisvy,NSTEP,NREC,DELTAT)
- call create_2D_image(vx(:,:,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx(:,:,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,1)
- call create_2D_image(vy(:,:,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy(:,:,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,2)
endif
@@ -1254,7 +1254,7 @@
write(20,*) 'set xlabel "Time (s)"'
write(20,*) 'set ylabel "Total energy"'
write(20,*)
- write(20,*) 'set output "collino3D_total_energy_semilog.eps"'
+ write(20,*) 'set output "CPML3D_total_energy_semilog.eps"'
write(20,*) 'set logscale y'
write(20,*) 'plot "energy.dat" t ''Total energy'' w l 1'
write(20,*) 'pause -1 "Hit any key..."'
@@ -1358,7 +1358,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -1391,6 +1391,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -1495,7 +1496,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_viscoelastic_MPI.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_viscoelastic_MPI.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_CPML_3D_viscoelastic_MPI.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -66,17 +66,6 @@
! number = {1},
! doi = {10.1111/j.1365-246X.2009.04278.x}}
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdela\^aziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -98,6 +87,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
@@ -212,8 +212,8 @@
double precision, parameter :: NPOWER = 2.d0
double precision, parameter :: K_MAX_PML = 7.d0 ! from Gedney page 8.11
-! double precision, parameter :: ALPHA_MAX_PML = 0.d0 ! from festa and Vilotte
- double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from festa and Vilotte
+! double precision, parameter :: ALPHA_MAX_PML = 0.d0 ! from Festa and Vilotte
+ double precision, parameter :: ALPHA_MAX_PML = 2.d0*PI*(f0/2.d0) ! from Festa and Vilotte
! arrays for the memory variables
! could declare these arrays in PML only to save a lot of memory, but proof of concept only here
@@ -259,9 +259,9 @@
double precision :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz,div
! 1D arrays for the damping profiles
- double precision, dimension(1:NX) :: d_x,K_x,alpha_prime_x,a_x,b_x,d_x_half,K_x_half,alpha_prime_x_half,a_x_half,b_x_half
- double precision, dimension(1:NY) :: d_y,K_y,alpha_prime_y,a_y,b_y,d_y_half,K_y_half,alpha_prime_y_half,a_y_half,b_y_half
- double precision, dimension(1:NZ) :: d_z,K_z,alpha_prime_z,a_z,b_z,d_z_half,K_z_half,alpha_prime_z_half,a_z_half,b_z_half
+ double precision, dimension(1:NX) :: d_x,K_x,alpha_x,a_x,b_x,d_x_half,K_x_half,alpha_x_half,a_x_half,b_x_half
+ double precision, dimension(1:NY) :: d_y,K_y,alpha_y,a_y,b_y,d_y_half,K_y_half,alpha_y_half,a_y_half,b_y_half
+ double precision, dimension(1:NZ) :: d_z,K_z,alpha_z,a_z,b_z,d_z_half,K_z_half,alpha_z_half,a_z_half,b_z_half
! PML
double precision thickness_PML_x,thickness_PML_y,thickness_PML_z
@@ -476,13 +476,13 @@
thickness_PML_y = NPOINTS_PML * DELTAY
thickness_PML_z = NPOINTS_PML * DELTAZ
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.0001d0
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0_x = - (NPOWER + 1) * cp *dsqrt(taumax)* log(Rcoef) / (2.d0 * thickness_PML_x)
d0_y = - (NPOWER + 1) * cp *dsqrt(taumax)* log(Rcoef) / (2.d0 * thickness_PML_y)
d0_z = - (NPOWER + 1) * cp *dsqrt(taumax)* log(Rcoef) / (2.d0 * thickness_PML_z)
@@ -499,8 +499,8 @@
d_x_half(:) = ZERO
K_x(:) = 1.d0
K_x_half(:) = 1.d0
- alpha_prime_x(:) = ZERO
- alpha_prime_x_half(:) = ZERO
+ alpha_x(:) = ZERO
+ alpha_x_half(:) = ZERO
a_x(:) = ZERO
a_x_half(:) = ZERO
@@ -508,8 +508,8 @@
d_y_half(:) = ZERO
K_y(:) = 1.d0
K_y_half(:) = 1.d0
- alpha_prime_y(:) = ZERO
- alpha_prime_y_half(:) = ZERO
+ alpha_y(:) = ZERO
+ alpha_y_half(:) = ZERO
a_y(:) = ZERO
a_y_half(:) = ZERO
@@ -517,8 +517,8 @@
d_z_half(:) = ZERO
K_z(:) = 1.d0
K_z_half(:) = 1.d0
- alpha_prime_z(:) = ZERO
- alpha_prime_z_half(:) = ZERO
+ alpha_z(:) = ZERO
+ alpha_z_half(:) = ZERO
a_z(:) = ZERO
a_z_half(:) = ZERO
@@ -543,7 +543,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -553,7 +553,7 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -568,7 +568,7 @@
d_x(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -578,22 +578,22 @@
d_x_half(i) = d0_x * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_x_half(i) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_x_half(i) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
! just in case, for -5 at the end
- if(alpha_prime_x(i) < ZERO) alpha_prime_x(i) = ZERO
- if(alpha_prime_x_half(i) < ZERO) alpha_prime_x_half(i) = ZERO
+ if(alpha_x(i) < ZERO) alpha_x(i) = ZERO
+ if(alpha_x_half(i) < ZERO) alpha_x_half(i) = ZERO
- b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_prime_x(i)) * DELTAT)
- b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_prime_x_half(i)) * DELTAT)
+ b_x(i) = exp(- (d_x(i) / K_x(i) + alpha_x(i)) * DELTAT)
+ b_x_half(i) = exp(- (d_x_half(i) / K_x_half(i) + alpha_x_half(i)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_prime_x(i)))
+ if(abs(d_x(i)) > 1.d-6) a_x(i) = d_x(i) * (b_x(i) - 1.d0) / (K_x(i) * (d_x(i) + K_x(i) * alpha_x(i)))
if(abs(d_x_half(i)) > 1.d-6) a_x_half(i) = d_x_half(i) * &
- (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_prime_x_half(i)))
+ (b_x_half(i) - 1.d0) / (K_x_half(i) * (d_x_half(i) + K_x_half(i) * alpha_x_half(i)))
enddo
@@ -618,7 +618,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -628,7 +628,7 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -643,7 +643,7 @@
d_y(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -653,18 +653,18 @@
d_y_half(j) = d0_y * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_y_half(j) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_y_half(j) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_prime_y(j)) * DELTAT)
- b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_prime_y_half(j)) * DELTAT)
+ b_y(j) = exp(- (d_y(j) / K_y(j) + alpha_y(j)) * DELTAT)
+ b_y_half(j) = exp(- (d_y_half(j) / K_y_half(j) + alpha_y_half(j)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_prime_y(j)))
+ if(abs(d_y(j)) > 1.d-6) a_y(j) = d_y(j) * (b_y(j) - 1.d0) / (K_y(j) * (d_y(j) + K_y(j) * alpha_y(j)))
if(abs(d_y_half(j)) > 1.d-6) a_y_half(j) = d_y_half(j) * &
- (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_prime_y_half(j)))
+ (b_y_half(j) - 1.d0) / (K_y_half(j) * (d_y_half(j) + K_y_half(j) * alpha_y_half(j)))
enddo
@@ -689,7 +689,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -699,7 +699,7 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
@@ -714,7 +714,7 @@
d_z(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
! define damping profile at half the grid points
@@ -724,18 +724,18 @@
d_z_half(k) = d0_z * abscissa_normalized**NPOWER
! this taken from Gedney page 8.2
K_z_half(k) = 1.d0 + (K_MAX_PML - 1.d0) * abscissa_normalized**NPOWER
- alpha_prime_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
+ alpha_z_half(k) = ALPHA_MAX_PML * (1.d0 - abscissa_normalized)
endif
endif
- b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_prime_z(k)) * DELTAT)
- b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_prime_z_half(k)) * DELTAT)
+ b_z(k) = exp(- (d_z(k) / K_z(k) + alpha_z(k)) * DELTAT)
+ b_z_half(k) = exp(- (d_z_half(k) / K_z_half(k) + alpha_z_half(k)) * DELTAT)
! this to avoid division by zero outside the PML
- if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_prime_z(k)))
+ if(abs(d_z(k)) > 1.d-6) a_z(k) = d_z(k) * (b_z(k) - 1.d0) / (K_z(k) * (d_z(k) + K_z(k) * alpha_z(k)))
if(abs(d_z_half(k)) > 1.d-6) a_z_half(k) = d_z_half(k) * &
- (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_prime_z_half(k)))
+ (b_z_half(k) - 1.d0) / (K_z_half(k) * (d_z_half(k) + K_z_half(k) * alpha_z_half(k)))
enddo
@@ -1436,9 +1436,9 @@
print *
call write_seismograms(sisvx,sisvy,NSTEP,NREC,DELTAT,t0)
- call create_2D_image(vx(1:NX,1:NY,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx(1:NX,1:NY,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,1,max_amplitudeVx)
- call create_2D_image(vy(1:NX,1:NY,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy(1:NX,1:NY,NZ_LOCAL),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,2,max_amplitudeVy)
endif
@@ -1569,7 +1569,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number,max_amplitude)
implicit none
@@ -1603,6 +1603,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
! write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -1707,7 +1708,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_anisotropic_fourth.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_anisotropic_fourth.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_anisotropic_fourth.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -50,20 +50,6 @@
implicit none
!
-! If you use this code for your own research, please cite some (or all) of these articles:
-!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-
-!
! PML implemented in the two directions (x and y directions).
!
! Version 1.0 July, 2010
@@ -107,7 +93,6 @@
integer, parameter :: NX = 401
integer, parameter :: NY = 401
-
! size of a grid cell
double precision, parameter :: h = 5.d0
@@ -244,7 +229,7 @@
print *,'Velocity of qSV along horizontal axis . . =',sqrt(c33/rho)
print *
-! from Becache et al., INRIA report, equation 7 page 5
+! from Becache et al., INRIA report, equation 7 page 5 http://hal.inria.fr/docs/00/07/22/83/PDF/RR-4304.pdf
if(c11*c22 - c12*c12 <= 0.d0) stop 'problem in definition of orthotropic material'
! check intrinsic mathematical stability of PML model for an anisotropic material
@@ -278,10 +263,10 @@
! thickness of the layer in meters
delta = NPOINTS_PML * h
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0 = 3.d0 * quasi_cp_max * log(1.d0/Rcoef) / (2.d0 * delta)
print *,'d0 = ',d0
@@ -625,11 +610,11 @@
if(velocnorm > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
image_data_2D = vx_1 + vx_2
- call create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,.true.,.true.,.true.,.true.,1)
image_data_2D = vy_1 + vy_2
- call create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,.true.,.true.,.true.,.true.,2)
endif
@@ -745,7 +730,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -778,6 +763,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -882,7 +868,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_isotropic.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_isotropic.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_2D_isotropic.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -46,7 +46,7 @@
! Application of the PML Absorbing Layer Model to the Linear
! Elastodynamic Problem in Anisotropic Heteregeneous Media
! INRIA Research Report RR-3471, August 1998
-! http://www.inria.fr/publications
+! http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
!
! and
!
@@ -207,10 +207,10 @@
! thickness of the layer in meters
delta = NPOINTS_PML * h
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0 = 3.d0 * cp * log(1.d0/Rcoef) / (2.d0 * delta)
print *,'d0 = ',d0
@@ -542,11 +542,11 @@
if(velocnorm > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
image_data_2D = vx_1 + vx_2
- call create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,.true.,.true.,.true.,.true.,1)
image_data_2D = vy_1 + vy_2
- call create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,.true.,.true.,.true.,.true.,2)
endif
@@ -662,7 +662,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -695,6 +695,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -799,7 +800,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_3D_isotropic_OpenMP.f90
===================================================================
--- seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_3D_isotropic_OpenMP.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/tags/v1.1.3/seismic_PML_Collino_3D_isotropic_OpenMP.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -47,7 +47,7 @@
! Application of the PML Absorbing Layer Model to the Linear
! Elastodynamic Problem in Anisotropic Heteregeneous Media
! INRIA Research Report RR-3471, August 1998
-! http://www.inria.fr/publications
+! http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
!
! and
!
@@ -208,10 +208,10 @@
! thickness of the layer in meters
delta = NPOINTS_PML * h
-! reflection coefficient (INRIA report section 6.1)
+! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
-! compute d0 from INRIA report section 6.1
+! compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
d0 = 3.d0 * cp * log(1.d0/Rcoef) / (2.d0 * delta)
print *,'d0 = ',d0
@@ -862,10 +862,10 @@
print *
call write_seismograms(sisvx,sisvy,NSTEP,NREC,DELTAT)
- call create_2D_image(vx_1(:,:,NZ/2) + vx_2(:,:,NZ/2) + vx_3(:,:,NZ/2),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vx_1(:,:,NZ/2) + vx_2(:,:,NZ/2) + vx_3(:,:,NZ/2),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,.true.,.true.,.true.,.true.,1)
- call create_2D_image(vy_1(:,:,NZ/2) + vy_2(:,:,NZ/2) +vy_3(:,:,NZ/2),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ call create_color_image(vy_1(:,:,NZ/2) + vy_2(:,:,NZ/2) +vy_3(:,:,NZ/2),NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,.true.,.true.,.true.,.true.,2)
endif
@@ -989,7 +989,7 @@
!---- the image is created in PNM format and then converted to GIF
!----
- subroutine create_2D_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
+ subroutine create_color_image(image_data_2D,NX,NY,it,ISOURCE,JSOURCE,ix_rec,iy_rec,nrec, &
NPOINTS_PML,USE_PML_XMIN,USE_PML_XMAX,USE_PML_YMIN,USE_PML_YMAX,field_number)
implicit none
@@ -1022,6 +1022,7 @@
double precision :: normalized_value,max_amplitude
! open image file and create system command to convert image to more convenient format
+! use the "convert" command from ImageMagick http://www.imagemagick.org
if(field_number == 1) then
write(file_name,"('image',i6.6,'_Vx.pnm')") it
write(system_command,"('convert image',i6.6,'_Vx.pnm image',i6.6,'_Vx.gif ; rm image',i6.6,'_Vx.pnm')") it,it,it
@@ -1126,7 +1127,7 @@
! call the system to convert image to GIF (can be commented out if "call system" is missing in your compiler)
! call system(system_command)
- end subroutine create_2D_image
+ end subroutine create_color_image
!
! CeCILL FREE SOFTWARE LICENSE AGREEMENT
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_anisotropic.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -71,17 +71,6 @@
! If you use this code for your own research, please cite some (or all) of these
! articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoGe08,
! author = {Roland Martin and Dimitri Komatitsch and Stephen D. Gedney},
! title = {A variational formulation of a stabilized unsplit convolutional perfectly
@@ -114,6 +103,17 @@
! number = {1},
! doi = {10.1111/j.1365-246X.2009.04278.x}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! If you use the anisotropic implementation, please cite this article,
! in which the anisotropic parameters are described, as well:
!
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_isotropic_fourth_order.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -77,17 +77,6 @@
!
! If you use this code for your own research, please cite some (or all) of these articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdelaaziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -120,6 +109,17 @@
! pages = {334-339},
! doi = {10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! To display the 2D results as color images, use:
!
! " display image*.gif " or " gimp image*.gif "
Modified: seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/trunk/seismic_CPML_2D_poroelastic_fourth_order.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -82,17 +82,6 @@
! number = {4},
! doi = {10.1190/1.2939484}}
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKo09,
! author = {Roland Martin and Dimitri Komatitsch},
! title = {An unsplit convolutional perfectly matched layer technique improved
@@ -114,6 +103,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
Modified: seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/trunk/seismic_CPML_3D_isotropic_MPI_OpenMP.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -53,17 +53,6 @@
! If you use this code for your own research, please cite some (or all) of these
! articles:
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdela\^aziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -96,6 +85,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
Modified: seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/trunk/seismic_CPML_3D_viscoelastic_MPI.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -66,17 +66,6 @@
! number = {1},
! doi = {10.1111/j.1365-246X.2009.04278.x}}
!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-!
! @ARTICLE{MaKoEz08,
! author = {Roland Martin and Dimitri Komatitsch and Abdela\^aziz Ezziani},
! title = {An unsplit convolutional perfectly matched layer improved at grazing
@@ -98,6 +87,17 @@
! pages = {274-304},
! number = {3}}
!
+! @ARTICLE{KoMa07,
+! author = {Dimitri Komatitsch and Roland Martin},
+! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
+! at grazing incidence for the seismic wave equation},
+! journal = {Geophysics},
+! year = {2007},
+! volume = {72},
+! number = {5},
+! pages = {SM155-SM167},
+! doi = {10.1190/1.2757586}}
+!
! The original CPML technique for Maxwell's equations is described in:
!
! @ARTICLE{RoGe00,
Modified: seismo/3D/CPML/trunk/seismic_PML_Collino_2D_anisotropic_fourth.f90
===================================================================
--- seismo/3D/CPML/trunk/seismic_PML_Collino_2D_anisotropic_fourth.f90 2011-05-29 00:10:34 UTC (rev 18484)
+++ seismo/3D/CPML/trunk/seismic_PML_Collino_2D_anisotropic_fourth.f90 2011-05-29 00:17:47 UTC (rev 18485)
@@ -50,20 +50,6 @@
implicit none
!
-! If you use this code for your own research, please cite some (or all) of these articles:
-!
-! @ARTICLE{KoMa07,
-! author = {Dimitri Komatitsch and Roland Martin},
-! title = {An unsplit convolutional {P}erfectly {M}atched {L}ayer improved
-! at grazing incidence for the seismic wave equation},
-! journal = {Geophysics},
-! year = {2007},
-! volume = {72},
-! number = {5},
-! pages = {SM155-SM167},
-! doi = {10.1190/1.2757586}}
-
-!
! PML implemented in the two directions (x and y directions).
!
! Version 1.0 July, 2010
@@ -107,7 +93,6 @@
integer, parameter :: NX = 401
integer, parameter :: NY = 401
-
! size of a grid cell
double precision, parameter :: h = 5.d0
More information about the CIG-COMMITS
mailing list