[cig-commits] r15578 - in seismo/2D/SPECFEM2D/trunk: . DATA
cmorency at geodynamics.org
cmorency at geodynamics.org
Sun Aug 23 09:28:37 PDT 2009
Author: cmorency
Date: 2009-08-23 09:28:36 -0700 (Sun, 23 Aug 2009)
New Revision: 15578
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90
seismo/2D/SPECFEM2D/trunk/gmat01.f90
seismo/2D/SPECFEM2D/trunk/plotpost.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt
Log:
Coupled elastic/poro fixed for unstructured meshes. Cleaned unused variables.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION 2009-08-23 16:28:36 UTC (rev 15578)
@@ -1,8 +1,8 @@
#source 1
-source_surf = .true. # source inside the medium or at the surface
-xs = 2000 #4000. #2000. # source location x in meters
-zs = 00 #2000. #00. # source location z in meters
-source_type = 2 # elastic force or acoustic pressure = 1 or moment tensor = 2
+source_surf = .false. # source inside the medium or at the surface
+xs = 0. # source location x in meters
+zs = -100. # source location z in meters
+source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 5.0 # dominant source frequency (Hz) if not Dirac or Heaviside
t0 = 0.0 # time shift when multi sources (if one source, must be zero)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-08-23 16:28:36 UTC (rev 15578)
@@ -1,26 +1,26 @@
# title of job, and file that contains interface data
-title = 2D SEG Tests
-interfacesfile = interfaces.dat
+title = Test for M2 UPPA
+interfacesfile = interfaces_M2_UPPA_curved.dat
# data concerning mesh, when generated using third-party app (more info in README)
read_external_mesh = .true.
-mesh_file = ./DATA/XSEG_small/mesh_file # file containing the mesh
-nodes_coords_file = ./DATA/XSEG_small/nodes_coords_file # file containing the nodes coordinates
-materials_file = ./DATA/XSEG_small/materials_file # file containing the material number for each element
-free_surface_file = ./DATA/XSEG_small/free_surface_file # file containing the free surface
-absorbing_surface_file = ./DATA/XSEG_small/absorbing_surface_file # file containing the absorbing surface
+mesh_file = ./DATA/unstructured_fluide_solide_test/mesh # file containing the mesh
+nodes_coords_file = ./DATA/unstructured_fluide_solide_test/nodes_coords # file containing the nodes coordinates
+materials_file = ./DATA/unstructured_fluide_solide_test/mat # file containing the material number for each element
+free_surface_file = ./DATA/unstructured_fluide_solide_test/surface_free # file containing the free surface
+absorbing_surface_file = ./DATA/unstructured_fluide_solide_test/surface_abs # file containing the absorbing surface
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
partitioning_method = 1 # ascending order = 1, Metis = 2, Scotch = 3
partitioning_strategy = 01110 #b{strat=m{asc=g,low=g,rat=1.0,type=h,vert=4}} #01110 # options concerning partitioning strategy.
-
+
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
-xmax = 16000.d0 # abscissa of right side of the model
-nx = 300 # number of elements along X
+xmax = 4000.d0 # abscissa of right side of the model
+nx = 80 # number of elements along X
ngnod = 4 # number of control nodes per element (4 or 9)
initialfield = .false. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
@@ -32,22 +32,22 @@
freq0 = 10 # frequency for viscous attenuation
# absorbing boundaries parameters
-absorbing_conditions = .true. # absorbing boundary active or not
-absorbbottom = .true.
+absorbing_conditions = .true. # absorbing boundary active or not
+absorbbottom = .true.
absorbright = .true.
absorbtop = .false.
absorbleft = .true.
# time step parameters
-nt = 30000 # total number of time steps
-deltat = 4d-4 # duration of a time step
+nt = 18000 # total number of time steps
+deltat = 0.25d-3 # duration of a time step
isolver = 1 # type of simulation 1=forward 2=adjoint + kernels
# source parameters
NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
-# constants for attenuation
+# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
## DK DK Qp and Qs can now vary in each spectral element and are therefore given in the same list
## DK DK list as rho, Vp and Vs at the end of this file
@@ -56,7 +56,7 @@
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
# receiver line parameters for seismograms
-seismotype = 2 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl 6=potential
+seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
save_forward = .false. # save the last frame, needed for adjoint simulation
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
@@ -66,13 +66,13 @@
# first receiver line
nrec = 100 # number of receivers
xdeb = 0. # first receiver x in meters
-zdeb = 0. # first receiver z in meters
-xfin = 15000. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 0. # last receiver z in meters (ignored if onlyone receiver)
+zdeb = -30. # first receiver z in meters
+xfin = 13000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = -30. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf = .false. # receivers inside the medium or at the surface
# display parameters
-NTSTEP_BETWEEN_OUTPUT_INFO = 500 # display frequency in time steps
+NTSTEP_BETWEEN_OUTPUT_INFO = 2000 # display frequency in time steps
output_postscript_snapshot = .false. # output Postscript snapshot of the results
output_color_image = .true. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
@@ -88,34 +88,19 @@
outputgrid = .false. # save the grid in a text file or not
# velocity and density models
-nbmodels = 22 # nb of different models
+nbmodels = 4 # nb of different models
# define models as I: (model_number,1,rho,Vp,Vs,0,0,Qp,Qs) or II: (model_number,2,rho,c11,c13,c33,c44,Qp,Qs) or III: (model_number,3,rhos,rhof,phi,c,kxx,kxz,kzz,Ks,Kf,Kfr,etaf,mufr,Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, & poroelastic models simultaneously
-1 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-2 1 1925.8d0 2050.d0 1183.41d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-3 1 2312.d0 4500.d0 2597.58d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-4 1 1986.3d0 2200.d0 1270.05d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-5 1 1881.4d0 1950.d0 1125.67d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-6 1 1808.1d0 1800.d0 1039.03d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-7 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-8 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-9 1 1986.3d0 2200.d0 1270.06d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-10 1 1904.0d0 2000.d0 1154.55d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-11 1 1857.9d0 1900.d0 1096.80d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-12 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-13 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-14 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-15 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-16 1 2056.6d0 2400.d0 1385.52d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-17 1 2116.0d0 2600.d0 1501.10d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-18 1 2087.6d0 2500.d0 1443.35d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-19 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-20 1 1857.9d0 1900.d0 1096.80d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-21 1 2087.6d0 2500.d0 1443.35d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-#22 3 1986.3 1040.3 0.4 2.0 1d-11 0.0 1d-11 5.341d9 5.341d8 5.341d8 0.0d-4 3.204d9 10.d0
-22 3 2200.d0 786.3d0 0.4 2.0 1d-11 0.0 1d-11 5.341d9 2d9 3d9 0.0d-4 3.204d9 10.d0
-#22 1 1634.52d0 2272.89d0 1472.71d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+1 1 1000.d0 1500.d0 0.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+2 1 2500.d0 5000.d0 2500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+4 3 2200.d0 786.3d0 0.4 2.0 1d-11 0.0 1d-11 5.341d9 2d9 3d9 0.0d-4 3.204d9 10.d0
+#4 1 2200.d0 2200.d0 1343.375d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 1 # nb of regions and model number for each
-1 300 1 160 1
+nbregions = 5 # nb of regions and model number for each
+1 80 1 20 1
+1 80 21 40 4
+1 80 41 60 3
+60 80 21 40 4
+30 40 50 60 2
Modified: seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2009-08-23 16:28:36 UTC (rev 15578)
@@ -1,100 +1,100 @@
-S0001 AA 0.0000000 0.0000000 0.0 0.0
-S0002 AA 151.5151515 0.0000000 0.0 0.0
-S0003 AA 303.0303030 0.0000000 0.0 0.0
-S0004 AA 454.5454545 0.0000000 0.0 0.0
-S0005 AA 606.0606061 0.0000000 0.0 0.0
-S0006 AA 757.5757576 0.0000000 0.0 0.0
-S0007 AA 909.0909091 0.0000000 0.0 0.0
-S0008 AA 1060.6060606 0.0000000 0.0 0.0
-S0009 AA 1212.1212121 0.0000000 0.0 0.0
-S0010 AA 1363.6363636 0.0000000 0.0 0.0
-S0011 AA 1515.1515152 0.0000000 0.0 0.0
-S0012 AA 1666.6666667 0.0000000 0.0 0.0
-S0013 AA 1818.1818182 0.0000000 0.0 0.0
-S0014 AA 1969.6969697 0.0000000 0.0 0.0
-S0015 AA 2121.2121212 0.0000000 0.0 0.0
-S0016 AA 2272.7272727 0.0000000 0.0 0.0
-S0017 AA 2424.2424242 0.0000000 0.0 0.0
-S0018 AA 2575.7575758 0.0000000 0.0 0.0
-S0019 AA 2727.2727273 0.0000000 0.0 0.0
-S0020 AA 2878.7878788 0.0000000 0.0 0.0
-S0021 AA 3030.3030303 0.0000000 0.0 0.0
-S0022 AA 3181.8181818 0.0000000 0.0 0.0
-S0023 AA 3333.3333333 0.0000000 0.0 0.0
-S0024 AA 3484.8484848 0.0000000 0.0 0.0
-S0025 AA 3636.3636364 0.0000000 0.0 0.0
-S0026 AA 3787.8787879 0.0000000 0.0 0.0
-S0027 AA 3939.3939394 0.0000000 0.0 0.0
-S0028 AA 4090.9090909 0.0000000 0.0 0.0
-S0029 AA 4242.4242424 0.0000000 0.0 0.0
-S0030 AA 4393.9393939 0.0000000 0.0 0.0
-S0031 AA 4545.4545455 0.0000000 0.0 0.0
-S0032 AA 4696.9696970 0.0000000 0.0 0.0
-S0033 AA 4848.4848485 0.0000000 0.0 0.0
-S0034 AA 5000.0000000 0.0000000 0.0 0.0
-S0035 AA 5151.5151515 0.0000000 0.0 0.0
-S0036 AA 5303.0303030 0.0000000 0.0 0.0
-S0037 AA 5454.5454545 0.0000000 0.0 0.0
-S0038 AA 5606.0606061 0.0000000 0.0 0.0
-S0039 AA 5757.5757576 0.0000000 0.0 0.0
-S0040 AA 5909.0909091 0.0000000 0.0 0.0
-S0041 AA 6060.6060606 0.0000000 0.0 0.0
-S0042 AA 6212.1212121 0.0000000 0.0 0.0
-S0043 AA 6363.6363636 0.0000000 0.0 0.0
-S0044 AA 6515.1515152 0.0000000 0.0 0.0
-S0045 AA 6666.6666667 0.0000000 0.0 0.0
-S0046 AA 6818.1818182 0.0000000 0.0 0.0
-S0047 AA 6969.6969697 0.0000000 0.0 0.0
-S0048 AA 7121.2121212 0.0000000 0.0 0.0
-S0049 AA 7272.7272727 0.0000000 0.0 0.0
-S0050 AA 7424.2424242 0.0000000 0.0 0.0
-S0051 AA 7575.7575758 0.0000000 0.0 0.0
-S0052 AA 7727.2727273 0.0000000 0.0 0.0
-S0053 AA 7878.7878788 0.0000000 0.0 0.0
-S0054 AA 8030.3030303 0.0000000 0.0 0.0
-S0055 AA 8181.8181818 0.0000000 0.0 0.0
-S0056 AA 8333.3333333 0.0000000 0.0 0.0
-S0057 AA 8484.8484848 0.0000000 0.0 0.0
-S0058 AA 8636.3636364 0.0000000 0.0 0.0
-S0059 AA 8787.8787879 0.0000000 0.0 0.0
-S0060 AA 8939.3939394 0.0000000 0.0 0.0
-S0061 AA 9090.9090909 0.0000000 0.0 0.0
-S0062 AA 9242.4242424 0.0000000 0.0 0.0
-S0063 AA 9393.9393939 0.0000000 0.0 0.0
-S0064 AA 9545.4545455 0.0000000 0.0 0.0
-S0065 AA 9696.9696970 0.0000000 0.0 0.0
-S0066 AA 9848.4848485 0.0000000 0.0 0.0
-S0067 AA 10000.0000000 0.0000000 0.0 0.0
-S0068 AA 10151.5151515 0.0000000 0.0 0.0
-S0069 AA 10303.0303030 0.0000000 0.0 0.0
-S0070 AA 10454.5454545 0.0000000 0.0 0.0
-S0071 AA 10606.0606061 0.0000000 0.0 0.0
-S0072 AA 10757.5757576 0.0000000 0.0 0.0
-S0073 AA 10909.0909091 0.0000000 0.0 0.0
-S0074 AA 11060.6060606 0.0000000 0.0 0.0
-S0075 AA 11212.1212121 0.0000000 0.0 0.0
-S0076 AA 11363.6363636 0.0000000 0.0 0.0
-S0077 AA 11515.1515152 0.0000000 0.0 0.0
-S0078 AA 11666.6666667 0.0000000 0.0 0.0
-S0079 AA 11818.1818182 0.0000000 0.0 0.0
-S0080 AA 11969.6969697 0.0000000 0.0 0.0
-S0081 AA 12121.2121212 0.0000000 0.0 0.0
-S0082 AA 12272.7272727 0.0000000 0.0 0.0
-S0083 AA 12424.2424242 0.0000000 0.0 0.0
-S0084 AA 12575.7575758 0.0000000 0.0 0.0
-S0085 AA 12727.2727273 0.0000000 0.0 0.0
-S0086 AA 12878.7878788 0.0000000 0.0 0.0
-S0087 AA 13030.3030303 0.0000000 0.0 0.0
-S0088 AA 13181.8181818 0.0000000 0.0 0.0
-S0089 AA 13333.3333333 0.0000000 0.0 0.0
-S0090 AA 13484.8484848 0.0000000 0.0 0.0
-S0091 AA 13636.3636364 0.0000000 0.0 0.0
-S0092 AA 13787.8787879 0.0000000 0.0 0.0
-S0093 AA 13939.3939394 0.0000000 0.0 0.0
-S0094 AA 14090.9090909 0.0000000 0.0 0.0
-S0095 AA 14242.4242424 0.0000000 0.0 0.0
-S0096 AA 14393.9393939 0.0000000 0.0 0.0
-S0097 AA 14545.4545455 0.0000000 0.0 0.0
-S0098 AA 14696.9696970 0.0000000 0.0 0.0
-S0099 AA 14848.4848485 0.0000000 0.0 0.0
-S0100 AA 15000.0000000 0.0000000 0.0 0.0
+S0001 AA 0.0000000 -30.0000000 0.0 0.0
+S0002 AA 131.3131313 -30.0000000 0.0 0.0
+S0003 AA 262.6262626 -30.0000000 0.0 0.0
+S0004 AA 393.9393939 -30.0000000 0.0 0.0
+S0005 AA 525.2525253 -30.0000000 0.0 0.0
+S0006 AA 656.5656566 -30.0000000 0.0 0.0
+S0007 AA 787.8787879 -30.0000000 0.0 0.0
+S0008 AA 919.1919192 -30.0000000 0.0 0.0
+S0009 AA 1050.5050505 -30.0000000 0.0 0.0
+S0010 AA 1181.8181818 -30.0000000 0.0 0.0
+S0011 AA 1313.1313131 -30.0000000 0.0 0.0
+S0012 AA 1444.4444444 -30.0000000 0.0 0.0
+S0013 AA 1575.7575758 -30.0000000 0.0 0.0
+S0014 AA 1707.0707071 -30.0000000 0.0 0.0
+S0015 AA 1838.3838384 -30.0000000 0.0 0.0
+S0016 AA 1969.6969697 -30.0000000 0.0 0.0
+S0017 AA 2101.0101010 -30.0000000 0.0 0.0
+S0018 AA 2232.3232323 -30.0000000 0.0 0.0
+S0019 AA 2363.6363636 -30.0000000 0.0 0.0
+S0020 AA 2494.9494949 -30.0000000 0.0 0.0
+S0021 AA 2626.2626263 -30.0000000 0.0 0.0
+S0022 AA 2757.5757576 -30.0000000 0.0 0.0
+S0023 AA 2888.8888889 -30.0000000 0.0 0.0
+S0024 AA 3020.2020202 -30.0000000 0.0 0.0
+S0025 AA 3151.5151515 -30.0000000 0.0 0.0
+S0026 AA 3282.8282828 -30.0000000 0.0 0.0
+S0027 AA 3414.1414141 -30.0000000 0.0 0.0
+S0028 AA 3545.4545455 -30.0000000 0.0 0.0
+S0029 AA 3676.7676768 -30.0000000 0.0 0.0
+S0030 AA 3808.0808081 -30.0000000 0.0 0.0
+S0031 AA 3939.3939394 -30.0000000 0.0 0.0
+S0032 AA 4070.7070707 -30.0000000 0.0 0.0
+S0033 AA 4202.0202020 -30.0000000 0.0 0.0
+S0034 AA 4333.3333333 -30.0000000 0.0 0.0
+S0035 AA 4464.6464646 -30.0000000 0.0 0.0
+S0036 AA 4595.9595960 -30.0000000 0.0 0.0
+S0037 AA 4727.2727273 -30.0000000 0.0 0.0
+S0038 AA 4858.5858586 -30.0000000 0.0 0.0
+S0039 AA 4989.8989899 -30.0000000 0.0 0.0
+S0040 AA 5121.2121212 -30.0000000 0.0 0.0
+S0041 AA 5252.5252525 -30.0000000 0.0 0.0
+S0042 AA 5383.8383838 -30.0000000 0.0 0.0
+S0043 AA 5515.1515152 -30.0000000 0.0 0.0
+S0044 AA 5646.4646465 -30.0000000 0.0 0.0
+S0045 AA 5777.7777778 -30.0000000 0.0 0.0
+S0046 AA 5909.0909091 -30.0000000 0.0 0.0
+S0047 AA 6040.4040404 -30.0000000 0.0 0.0
+S0048 AA 6171.7171717 -30.0000000 0.0 0.0
+S0049 AA 6303.0303030 -30.0000000 0.0 0.0
+S0050 AA 6434.3434343 -30.0000000 0.0 0.0
+S0051 AA 6565.6565657 -30.0000000 0.0 0.0
+S0052 AA 6696.9696970 -30.0000000 0.0 0.0
+S0053 AA 6828.2828283 -30.0000000 0.0 0.0
+S0054 AA 6959.5959596 -30.0000000 0.0 0.0
+S0055 AA 7090.9090909 -30.0000000 0.0 0.0
+S0056 AA 7222.2222222 -30.0000000 0.0 0.0
+S0057 AA 7353.5353535 -30.0000000 0.0 0.0
+S0058 AA 7484.8484848 -30.0000000 0.0 0.0
+S0059 AA 7616.1616162 -30.0000000 0.0 0.0
+S0060 AA 7747.4747475 -30.0000000 0.0 0.0
+S0061 AA 7878.7878788 -30.0000000 0.0 0.0
+S0062 AA 8010.1010101 -30.0000000 0.0 0.0
+S0063 AA 8141.4141414 -30.0000000 0.0 0.0
+S0064 AA 8272.7272727 -30.0000000 0.0 0.0
+S0065 AA 8404.0404040 -30.0000000 0.0 0.0
+S0066 AA 8535.3535354 -30.0000000 0.0 0.0
+S0067 AA 8666.6666667 -30.0000000 0.0 0.0
+S0068 AA 8797.9797980 -30.0000000 0.0 0.0
+S0069 AA 8929.2929293 -30.0000000 0.0 0.0
+S0070 AA 9060.6060606 -30.0000000 0.0 0.0
+S0071 AA 9191.9191919 -30.0000000 0.0 0.0
+S0072 AA 9323.2323232 -30.0000000 0.0 0.0
+S0073 AA 9454.5454545 -30.0000000 0.0 0.0
+S0074 AA 9585.8585859 -30.0000000 0.0 0.0
+S0075 AA 9717.1717172 -30.0000000 0.0 0.0
+S0076 AA 9848.4848485 -30.0000000 0.0 0.0
+S0077 AA 9979.7979798 -30.0000000 0.0 0.0
+S0078 AA 10111.1111111 -30.0000000 0.0 0.0
+S0079 AA 10242.4242424 -30.0000000 0.0 0.0
+S0080 AA 10373.7373737 -30.0000000 0.0 0.0
+S0081 AA 10505.0505051 -30.0000000 0.0 0.0
+S0082 AA 10636.3636364 -30.0000000 0.0 0.0
+S0083 AA 10767.6767677 -30.0000000 0.0 0.0
+S0084 AA 10898.9898990 -30.0000000 0.0 0.0
+S0085 AA 11030.3030303 -30.0000000 0.0 0.0
+S0086 AA 11161.6161616 -30.0000000 0.0 0.0
+S0087 AA 11292.9292929 -30.0000000 0.0 0.0
+S0088 AA 11424.2424242 -30.0000000 0.0 0.0
+S0089 AA 11555.5555556 -30.0000000 0.0 0.0
+S0090 AA 11686.8686869 -30.0000000 0.0 0.0
+S0091 AA 11818.1818182 -30.0000000 0.0 0.0
+S0092 AA 11949.4949495 -30.0000000 0.0 0.0
+S0093 AA 12080.8080808 -30.0000000 0.0 0.0
+S0094 AA 12212.1212121 -30.0000000 0.0 0.0
+S0095 AA 12343.4343434 -30.0000000 0.0 0.0
+S0096 AA 12474.7474747 -30.0000000 0.0 0.0
+S0097 AA 12606.0606061 -30.0000000 0.0 0.0
+S0098 AA 12737.3737374 -30.0000000 0.0 0.0
+S0099 AA 12868.6868687 -30.0000000 0.0 0.0
+S0100 AA 13000.0000000 -30.0000000 0.0 0.0
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2009-08-23 16:28:36 UTC (rev 15578)
@@ -62,19 +62,20 @@
# NOTE FOR USERS OF IFORT 10.0 AND ABOVE :
# Use of option -heap-arrays <size> can be required, depending on the size of the simulation.
# Another workaround can be to increase your stack size (ulimit -s).
-#F90 = ifort
+F90 = ifort
#F90 = mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
-#CC = gcc
-#FLAGS_NOCHECK=-O3 -xP -vec-report0 -e95 -std95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
-#FLAGS_CHECK = $(FLAGS_NOCHECK)
+CC = gcc
+FLAGS_NOCHECK=-O3 -xP -vec-report0 -e95 -std95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
+FLAGS_CHECK = $(FLAGS_NOCHECK)
# GNU gfortran
-F90 = gfortran
+#F90 = gfortran
#F90 = mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
-CC = gcc
+#CC = gcc
##FLAGS_NOCHECK = -O3 -march=opteron -m64 -mfpmath=sse,387
-FLAGS_NOCHECK = -std=f95 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Wunused -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -mcmodel=medium
-FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
+#FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O3 -pedantic -pedantic-errors -Wunused -Waliasing -Wampersand -Wline-truncation -Wsurprising -Wunderflow -fno-trapping-math # -mcmodel=medium
+##FLAGS_NOCHECK = -std=f95 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Wunused -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -mcmodel=medium
+#FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
# IBM
#####F90 = xlf_r
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2009-08-23 16:28:36 UTC (rev 15578)
@@ -121,7 +121,7 @@
real(kind=CUSTOM_REAL) :: kappal_f,rhol_f
real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr,phil,tortl
real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,rhol_bar
- real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G,mul_C,lambdal_C,lambdalplus2mul_C,mul_M,lambdal_M,lambdalplus2mul_M
+ real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
kinetic_energy = ZERO
potential_energy = ZERO
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90 2009-08-23 16:28:36 UTC (rev 15578)
@@ -43,13 +43,13 @@
!
!========================================================================
- subroutine compute_forces_fluid(npoin,nspec,myrank,nelemabs,numat,iglob_source, &
+ subroutine compute_forces_fluid(npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,angleforce,deltatcube, &
+ initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,poroelastic,codeabs, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
- b_accelw_poroelastic,b_velocw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
+ b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
@@ -70,7 +70,7 @@
include "constants.h"
integer :: NSOURCE, i_source
- integer, dimension(NSOURCE) ::iglob_source,ispec_selected_source,source_type,is_proc_source
+ integer, dimension(NSOURCE) ::ispec_selected_source,source_type,is_proc_source
integer :: npoin,nspec,nelemabs,numat,it,NSTEP
integer :: nrec,isolver,myrank
integer, dimension(nrec) :: ispec_selected_rec,which_proc_receiver
@@ -84,7 +84,6 @@
logical :: save_forward
double precision ::deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
- double precision, dimension(NSOURCE) :: angleforce
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: kmato
@@ -96,8 +95,7 @@
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: accelw_poroelastic,velocw_poroelastic,displw_poroelastic,&
displs_poroelastic,velocs_poroelastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
- b_velocw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic
double precision, dimension(2,numat) :: density
double precision, dimension(3,numat) :: permeability
double precision, dimension(numat) :: porosity,tortuosity
@@ -124,8 +122,8 @@
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
! viscous attenuation
- double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous,b_rx_viscous
- double precision, dimension(NGLLX,NGLLZ,nspec) :: rz_viscous,b_rz_viscous
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: rz_viscous
double precision :: theta_e,theta_s
logical TURN_VISCATTENUATION_ON
double precision, dimension(3):: bl_unrelaxed,bl_relaxed
@@ -152,12 +150,8 @@
real(kind=CUSTOM_REAL) :: dwx_dxl,dwz_dxl,dwx_dzl,dwz_dzl
real(kind=CUSTOM_REAL) :: b_dux_dxi,b_dux_dgamma,b_duz_dxi,b_duz_dgamma
real(kind=CUSTOM_REAL) :: b_dux_dxl,b_duz_dxl,b_dux_dzl,b_duz_dzl
- real(kind=CUSTOM_REAL) :: dsxx,dsxz,dszz
- real(kind=CUSTOM_REAL) :: b_dsxx,b_dsxz,b_dszz
real(kind=CUSTOM_REAL) :: b_dwx_dxi,b_dwx_dgamma,b_dwz_dxi,b_dwz_dgamma
real(kind=CUSTOM_REAL) :: b_dwx_dxl,b_dwz_dxl,b_dwx_dzl,b_dwz_dzl
- real(kind=CUSTOM_REAL) :: dwxx,dwxz,dwzz
- real(kind=CUSTOM_REAL) :: b_dwxx,b_dwxz,b_dwzz
real(kind=CUSTOM_REAL) :: sigma_xx,sigma_xz,sigma_zz
real(kind=CUSTOM_REAL) :: sigmap
real(kind=CUSTOM_REAL) :: b_sigma_xx,b_sigma_xz,b_sigma_zz
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90 2009-08-23 16:28:36 UTC (rev 15578)
@@ -43,13 +43,13 @@
!
!========================================================================
- subroutine compute_forces_solid(npoin,nspec,myrank,nelemabs,numat,iglob_source, &
+ subroutine compute_forces_solid(npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,angleforce,deltatcube, &
+ initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,poroelastic,codeabs, &
accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displw_poroelastic,&
- b_accels_poroelastic,b_displs_poroelastic,b_velocw_poroelastic,b_displw_poroelastic,&
+ b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
@@ -70,7 +70,7 @@
include "constants.h"
integer :: NSOURCE, i_source
- integer, dimension(NSOURCE) :: iglob_source,ispec_selected_source,source_type,is_proc_source
+ integer, dimension(NSOURCE) :: ispec_selected_source,source_type,is_proc_source
integer :: npoin,nspec,nelemabs,numat,it,NSTEP
integer :: nrec,isolver,myrank
integer, dimension(nrec) :: ispec_selected_rec,which_proc_receiver
@@ -84,7 +84,6 @@
logical :: save_forward
double precision :: deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
- double precision, dimension(NSOURCE) :: angleforce
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: kmato
@@ -96,8 +95,7 @@
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: accels_poroelastic,velocs_poroelastic,displs_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: velocw_poroelastic,displw_poroelastic
- real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
- b_velocw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic
double precision, dimension(2,numat) :: density
double precision, dimension(3,numat) :: permeability
double precision, dimension(numat) :: porosity,tortuosity
@@ -124,8 +122,8 @@
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
! viscous attenuation (poroelastic media)
- double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous,b_rx_viscous
- double precision, dimension(NGLLX,NGLLZ,nspec) :: rz_viscous,b_rz_viscous
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: rz_viscous
double precision :: theta_e,theta_s
logical TURN_VISCATTENUATION_ON
double precision, dimension(3):: bl_unrelaxed,bl_relaxed
Modified: seismo/2D/SPECFEM2D/trunk/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gmat01.f90 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/gmat01.f90 2009-08-23 16:28:36 UTC (rev 15578)
@@ -62,7 +62,7 @@
integer in,n,indic
double precision young,poisson,cp,cs,mu,two_mu,lambda,Qp,Qs
double precision lambdaplus2mu_s,lambdaplus2mu_fr,kappa_s,kappa_f,kappa_fr
- double precision young_s,poisson_s,density(2),phi,tortuosity,permxx,permzz,permxz
+ double precision young_s,poisson_s,density(2),phi,tortuosity
double precision cpIsquare,cpIIsquare,cssquare,mu_s,mu_fr,eta_f,lambda_s,lambda_fr
double precision val1,val2,val3,val4,val5,val6
double precision val7,val8,val9,val10,val11,val12,val0
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.F90 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.F90 2009-08-23 16:28:36 UTC (rev 15578)
@@ -133,7 +133,7 @@
equivalence (postscript_line,ch1)
logical :: first
- double precision convert,x1,rlamda,rmu,denst,rKvol,cpIloc,xa,za,xb,zb
+ double precision convert,x1,cpIloc,xa,za,xb,zb
double precision z1,x2,z2,d,d1,d2,dummy,theta,thetaup,thetadown
double precision :: mul_s,kappal_s,rhol_s
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-08-23 16:28:36 UTC (rev 15578)
@@ -112,8 +112,9 @@
! journal={Geophys. J. Int.},
! doi={10.1111/j.1365-246X.2009.04332}}
!
-! version 6.0, Christina Morency, August 2009:
-! - support for poroelastic media (for regular meshes only for now).
+! version 6.0, Christina Morency and Yang Luo, August 2009:
+! - support for poroelastic media
+! - adjoint method for acoustic/elastic/poroelastic
!
! version 5.2, Dimitri Komatitsch, Nicolas Le Goff and Roland Martin, February 2008:
! - MPI implementation of the code based on domain decomposition
@@ -333,7 +334,7 @@
double precision :: permlxx,permlxz,permlzz,invpermlxx,invpermlxz,invpermlzz,detk
! adjoint
double precision, dimension(:), allocatable :: b_viscodampx,b_viscodampz
- integer reclen,reclen1,reclen2
+ integer reclen
! for fluid/solid coupling and edge detection
logical, dimension(:), allocatable :: elastic
@@ -354,7 +355,7 @@
integer, dimension(:), allocatable :: fluid_poro_acoustic_ispec,fluid_poro_acoustic_iedge, &
fluid_poro_poroelastic_ispec,fluid_poro_poroelastic_iedge
integer :: fluid_poro_acoustic_ispec_read, fluid_poro_poro_ispec_read
- integer :: num_fluid_poro_edges,num_fluid_poro_edges_alloc,iedge_poroelastic
+ integer :: num_fluid_poro_edges,iedge_poroelastic
logical :: coupled_acoustic_poro
double precision :: mul_G,lambdal_G,lambdalplus2mul_G
double precision :: dux_dxi,dux_dgamma,duz_dxi,duz_dgamma
@@ -370,17 +371,14 @@
integer, dimension(:), allocatable :: solid_poro_elastic_ispec,solid_poro_elastic_iedge, &
solid_poro_poroelastic_ispec,solid_poro_poroelastic_iedge
integer :: solid_poro_elastic_ispec_read, solid_poro_poro_ispec_read
- integer :: num_solid_poro_edges,num_solid_poro_edges_alloc,ispec_poroelastic,ii2,jj2
+ integer :: num_solid_poro_edges,ispec_poroelastic,ii2,jj2
logical :: coupled_elastic_poro
- double precision, dimension(:,:), allocatable :: displ,veloc
+ integer, dimension(:), allocatable :: icount
double precision :: sigma_xx,sigma_xz,sigma_zz,sigmap
double precision :: b_sigma_xx,b_sigma_xz,b_sigma_zz,b_sigmap
integer, dimension(:), allocatable :: ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,&
iend_top_poro,jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro
-! for detection elastic and acoustic valences
- integer, dimension(:), allocatable :: valence_elastic,valence_acoustic,valence_poroelastic
-
! for adjoint method
logical :: save_forward ! whether or not the last frame is saved to reconstruct the forward field
integer :: isolver ! 1 = forward wavefield, 2 = backward and adjoint wavefields and kernels
@@ -405,12 +403,11 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: phil_global,etal_f_global,rhol_s_global,rhol_f_global,rhol_bar_global, &
tortl_global,mulfr_global
real(kind=CUSTOM_REAL), dimension(:), allocatable :: permlxx_global,permlxz_global,permlzz_global
- character(len=150) :: adj_source_file,filename,filename2,filename3
+ character(len=150) :: adj_source_file
integer :: irec_local,nadj_rec_local
double precision :: xx,zz,rholb
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: adj_sourcearray
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: adj_sourcearrays
- double precision :: rhopmin,rhopmax,alphamin,alphamax,betamin,betamax
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_absorb_elastic_left,b_absorb_poro_s_left,b_absorb_poro_w_left
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_absorb_elastic_right,b_absorb_poro_s_right,b_absorb_poro_w_right
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_absorb_elastic_bottom,b_absorb_poro_s_bottom,b_absorb_poro_w_bottom
@@ -2142,11 +2139,9 @@
endif
if(any_poroelastic .and. any_elastic) then
- allocate(displ(2,npoin))
- allocate(veloc(2,npoin))
+ allocate(icount(npoin))
else
- allocate(displ(2,1))
- allocate(veloc(2,1))
+ allocate(icount(1))
endif
! potential, its first and second derivative, and inverse of the mass matrix for acoustic elements
@@ -4196,50 +4191,6 @@
endif
-! detecting poroelastic, elastic and acoustic global points valence
-
- if(coupled_acoustic_elastic .or. coupled_acoustic_poro .or. coupled_elastic_poro)then
-
- allocate(valence_elastic(npoin))
- allocate(valence_poroelastic(npoin))
- allocate(valence_acoustic(npoin))
-
- valence_elastic(:) = 0
- valence_poroelastic(:) = 0
- valence_acoustic(:) = 0
- do ispec = 1,nspec
- if(elastic(ispec)) then ! the element is elastic
- do k = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,k,ispec)
- valence_elastic(iglob) = valence_elastic(iglob) + 1
- enddo
- enddo
- elseif(poroelastic(ispec)) then ! the element is poroelastic
- do k = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,k,ispec)
- valence_poroelastic(iglob) = valence_poroelastic(iglob) + 1
- enddo
- enddo
- else ! the element is acoustic
- do k = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,k,ispec)
- valence_acoustic(iglob) = valence_acoustic(iglob) + 1
- enddo
- enddo
- endif
- enddo !do ispec
-
- else
-
- allocate(valence_elastic(1))
- allocate(valence_poroelastic(1))
- allocate(valence_acoustic(1))
-
- endif !(coupled_acoustic_elastic .or. coupled_acoustic_poro .or. coupled_elastic_poro)
-
#ifdef USE_MPI
if(OUTPUT_ENERGY) stop 'energy calculation only currently serial only, should add an MPI_REDUCE in parallel'
#endif
@@ -4820,21 +4771,11 @@
! compute dot product
displ_n = displ_x*nx + displ_z*nz
- if(valence_acoustic(iglob) == valence_elastic(iglob)) then
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + weight*displ_n
- else
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + &
- weight*displ_n*valence_acoustic(iglob)
- endif
if(isolver == 2) then
- if(valence_acoustic(iglob) == valence_elastic(iglob)) then
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) +&
weight*(b_displ_x*nx + b_displ_z*nz)
- else
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) + &
- weight*(b_displ_x*nx + b_displ_z*nz)*valence_acoustic(iglob)
- endif
endif !if(isolver == 2) then
enddo
@@ -4927,22 +4868,11 @@
! compute dot product [u_s + w]*n
displ_n = (displ_x + displw_x)*nx + (displ_z + displw_z)*nz
- if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + weight*displ_n
- else
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + &
- weight*displ_n*valence_acoustic(iglob)
- endif
if(isolver == 2) then
- if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) +&
weight*((b_displ_x + b_displw_x)*nx + (b_displ_z + b_displw_z)*nz)
- else
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) + &
- weight*((b_displ_x + b_displw_x)*nx + (b_displ_z + b_displw_z)*nz)* &
- valence_acoustic(iglob)
- endif
endif !if(isolver == 2) then
enddo
@@ -5050,86 +4980,6 @@
endif
-! ****************************************************************************************
-! If coupling elastic/poroelastic domain, average some arrays at the interface first
-! ****************************************************************************************
- if(coupled_elastic_poro) then
-
-! loop on all the coupling edges
- do inum = 1,num_solid_poro_edges
-
-! get the edge of the elastic element
- ispec_elastic = solid_poro_elastic_ispec(inum)
- iedge_elastic = solid_poro_elastic_iedge(inum)
-
-! get the corresponding edge of the poroelastic element
- ispec_poroelastic = solid_poro_poroelastic_ispec(inum)
- iedge_poroelastic = solid_poro_poroelastic_iedge(inum)
- phil = porosity(kmato(ispec_poroelastic))
-! implement 1D coupling along the edge
- do ipoin1D = 1,NGLLX
-
-! get point values for the poroelastic side, which matches our side in the inverse direction
- i = ivalue(ipoin1D,iedge_poroelastic)
- j = jvalue(ipoin1D,iedge_poroelastic)
- iglob = ibool(i,j,ispec_poroelastic)
-
-! get point values for the elastic side
- ii2 = ivalue_inverse(ipoin1D,iedge_elastic)
- jj2 = jvalue_inverse(ipoin1D,iedge_elastic)
- iglob2 = ibool(ii2,jj2,ispec_elastic)
-
- if(iglob /= iglob2) &
- call exit_MPI( 'error in solid/porous iglob detection')
-
- displ(1,iglob)=(2.d0*displs_poroelastic(1,iglob) + &
- displ_elastic(1,iglob2))/3.d0
- displ(2,iglob)=(2.d0*displs_poroelastic(2,iglob) + &
- displ_elastic(2,iglob2))/3.d0
-
- veloc(1,iglob)=(2.d0*velocs_poroelastic(1,iglob) +&
- veloc_elastic(1,iglob2))/3.d0
- veloc(2,iglob)=(2.d0*velocs_poroelastic(2,iglob) +&
- veloc_elastic(2,iglob2))/3.d0
-
- enddo
- enddo
-
-! loop on all the coupling edges
- do inum = 1,num_solid_poro_edges
-
-! get the corresponding edge of the poroelastic element
- ispec_poroelastic = solid_poro_poroelastic_ispec(inum)
- iedge_poroelastic = solid_poro_poroelastic_iedge(inum)
-! imnplement 1D coupling along the edge
- do ipoin1D = 1,NGLLX
-! get point values for the poroelastic side, which matches our side in the inverse direction
- i = ivalue(ipoin1D,iedge_poroelastic)
- j = jvalue(ipoin1D,iedge_poroelastic)
- iglob = ibool(i,j,ispec_poroelastic)
-
- displs_poroelastic(1,iglob)=displ(1,iglob)
- displs_poroelastic(2,iglob)=displ(2,iglob)
-
- displ_elastic(1,iglob)=displ(1,iglob)
- displ_elastic(2,iglob)=displ(2,iglob)
-
- velocs_poroelastic(1,iglob)=veloc(1,iglob)
- velocs_poroelastic(2,iglob)=veloc(2,iglob)
-
- veloc_elastic(1,iglob)=veloc(1,iglob)
- veloc_elastic(2,iglob)=veloc(2,iglob)
-
- displw_poroelastic(1,iglob)=ZERO
- displw_poroelastic(2,iglob)=ZERO
-
- velocw_poroelastic(1,iglob)=ZERO
- velocw_poroelastic(2,iglob)=ZERO
- enddo
- enddo
-
- endif
-
! *********************************************************
! ************* main solver for the elastic elements
! *********************************************************
@@ -5273,26 +5123,12 @@
weight = jacobian1D * wzgll(j)
endif
- if(valence_acoustic(iglob) == valence_elastic(iglob)) then
accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure
accel_elastic(2,iglob) = accel_elastic(2,iglob) + weight*nz*pressure
- else
- accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure*&
- valence_elastic(iglob)
- accel_elastic(2,iglob) = accel_elastic(2,iglob) + weight*nz*pressure*&
- valence_elastic(iglob)
- endif
if(isolver == 2) then
- if(valence_acoustic(iglob) == valence_elastic(iglob)) then
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) + weight*nz*b_pressure
- else
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure*&
- valence_elastic(iglob)
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) + weight*nz*b_pressure*&
- valence_elastic(iglob)
- endif
endif !if(isolver == 2) then
enddo
@@ -5562,34 +5398,18 @@
weight = jacobian1D * wzgll(j)
endif
- if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
accel_elastic(1,iglob) = accel_elastic(1,iglob) - weight* &
- (sigma_xx*nx + sigma_xz*nz +phil*sigmap*nx)/3.d0
+ (sigma_xx*nx + sigma_xz*nz +sigmap*nx)/3.d0
accel_elastic(2,iglob) = accel_elastic(2,iglob) - weight* &
- (sigma_xz*nx + sigma_zz*nz +phil*sigmap*nz)/3.d0
- else
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - weight* &
- (sigma_xx*nx + sigma_xz*nz +phil*sigmap*nx)/3.d0*valence_elastic(iglob)
+ (sigma_xz*nx + sigma_zz*nz +sigmap*nz)/3.d0
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - weight* &
- (sigma_xz*nx + sigma_zz*nz +phil*sigmap*nz)/3.d0*valence_elastic(iglob)
- endif
-
if(isolver == 2) then
- if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - weight* &
- (b_sigma_xx*nx + b_sigma_xz*nz +phil*b_sigmap*nx)/3.d0
+ (b_sigma_xx*nx + b_sigma_xz*nz +b_sigmap*nx)/3.d0
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - weight* &
- (b_sigma_xz*nx + b_sigma_zz*nz +phil*b_sigmap*nz)/3.d0
- else
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - weight* &
- (b_sigma_xx*nx + b_sigma_xz*nz +phil*b_sigmap*nx)/3.d0*valence_elastic(iglob)
-
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - weight* &
- (b_sigma_xz*nx + b_sigma_zz*nz +phil*b_sigmap*nz)/3.d0*valence_elastic(iglob)
- endif
+ (b_sigma_xz*nx + b_sigma_zz*nz +b_sigmap*nz)/3.d0
endif !if(isolver == 2) then
enddo
@@ -5689,13 +5509,13 @@
b_viscodampz(:) = ZERO
endif
- call compute_forces_solid(npoin,nspec,myrank,nelemabs,numat,iglob_source, &
+ call compute_forces_solid(npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,angleforce,deltatcube, &
+ initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,poroelastic,codeabs, &
accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displw_poroelastic,&
- b_accels_poroelastic,b_displs_poroelastic,b_velocw_poroelastic,b_displw_poroelastic,&
+ b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
@@ -5712,13 +5532,13 @@
- call compute_forces_fluid(npoin,nspec,myrank,nelemabs,numat,iglob_source, &
+ call compute_forces_fluid(npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,angleforce,deltatcube, &
+ initialfield,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,poroelastic,codeabs, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
- b_accelw_poroelastic,b_velocw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
+ b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
@@ -5872,7 +5692,6 @@
weight = jacobian1D * wzgll(j)
endif
- if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
! contribution to the solid phase
accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-phil/tortl)
accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-phil/tortl)
@@ -5880,22 +5699,8 @@
! contribution to the fluid phase
accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
- else
-! contribution to the solid phase
- accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-phil/tortl)*&
- valence_poroelastic(iglob)
- accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-phil/tortl)*&
- valence_poroelastic(iglob)
-! contribution to the fluid phase
- accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
- valence_poroelastic(iglob)
- accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
- valence_poroelastic(iglob)
- endif
-
if(isolver == 2) then
- if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
! contribution to the solid phase
b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-phil/tortl)
b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-phil/tortl)
@@ -5903,19 +5708,6 @@
! contribution to the fluid phase
b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
- else
-! contribution to the solid phase
- b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-phil/tortl)*&
- valence_poroelastic(iglob)
- b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-phil/tortl)*&
- valence_poroelastic(iglob)
-
-! contribution to the fluid phase
- b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
- valence_poroelastic(iglob)
- b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
- valence_poroelastic(iglob)
- endif
endif !if(isolver == 2) then
enddo ! do ipoin1D = 1,NGLLX
@@ -6188,54 +5980,26 @@
weight = jacobian1D * wzgll(j)
endif
- if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
! contribution to the solid phase
accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + &
- weight*((sigma_xx+phil*sigmap)*nx + sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)
+ weight*(sigma_xx*nx + sigma_xz*nz + sigmap*nx)/3.d0*(1.d0 -phil/tortl)
accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + &
- weight*(sigma_xz*nx + (sigma_zz+phil*sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)
+ weight*(sigma_xz*nx + sigma_zz*nz + sigmap*nz)/3.d0*(1.d0 -phil/tortl)
! contribution to the fluid phase
! w = 0
- else
-! contribution to the solid phase
- accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + &
- weight*((sigma_xx+phil*sigmap)*nx + sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)&
- *valence_poroelastic(iglob)
- accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + &
- weight*(sigma_xz*nx + (sigma_zz+phil*sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)&
- *valence_poroelastic(iglob)
-
-! contribution to the fluid phase
-! w = 0
- endif
-
if(isolver == 2) then
- if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
! contribution to the solid phase
b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + &
- weight*((b_sigma_xx+phil*b_sigmap)*nx + b_sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)
+ weight*(b_sigma_xx*nx + b_sigma_xz*nz + b_sigmap*nx)/3.d0*(1.d0 -phil/tortl)
b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + &
- weight*(b_sigma_xz*nx + (b_sigma_zz+phil*b_sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)
+ weight*(b_sigma_xz*nx + b_sigma_zz*nz + b_sigmap*nz)/3.d0*(1.d0 -phil/tortl)
! contribution to the fluid phase
! w = 0
- else
-! contribution to the solid phase
- b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + &
- weight*((b_sigma_xx+phil*b_sigmap)*nx + b_sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)&
- *valence_poroelastic(iglob)
-
- b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + &
- weight*(b_sigma_xz*nx + (b_sigma_zz+phil*b_sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)&
- *valence_poroelastic(iglob)
-
-! contribution to the fluid phase
-! w = 0
- endif
endif !if(isolver == 2) then
enddo
@@ -6357,6 +6121,61 @@
enddo
endif
+!assembling the displacements on the elastic-poro boundaries
+ if(coupled_elastic_poro) then
+ icount(:)=ZERO
+
+! loop on all the coupling edges
+ do inum = 1,num_solid_poro_edges
+! get the edge of the elastic element
+ ispec_elastic = solid_poro_elastic_ispec(inum)
+ iedge_elastic = solid_poro_elastic_iedge(inum)
+! get the corresponding edge of the poroelastic element
+ ispec_poroelastic = solid_poro_poroelastic_ispec(inum)
+ iedge_poroelastic = solid_poro_poroelastic_iedge(inum)
+
+ do ipoin1D = 1,NGLLX
+! recovering original velocities and accelerations on boundaries (elastic side)
+ i = ivalue_inverse(ipoin1D,iedge_elastic)
+ j = jvalue_inverse(ipoin1D,iedge_elastic)
+ iglob = ibool(i,j,ispec_elastic)
+ icount(iglob) = icount(iglob) + 1
+
+ if(icount(iglob) ==1)then
+ veloc_elastic(1,iglob) = veloc_elastic(1,iglob) - deltatover2*accel_elastic(1,iglob)
+ veloc_elastic(2,iglob) = veloc_elastic(2,iglob) - deltatover2*accel_elastic(2,iglob)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) / rmass_inverse_elastic(iglob)
+! recovering original velocities and accelerations on boundaries (poro side)
+ velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) - deltatover2*accels_poroelastic(1,iglob)
+ velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) - deltatover2*accels_poroelastic(2,iglob)
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
+! assembling accelerations
+ accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
+ ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ accel_elastic(2,iglob) = ( accel_elastic(2,iglob) + accels_poroelastic(2,iglob) ) / &
+ ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+ accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
+ accels_poroelastic(2,iglob) = accel_elastic(2,iglob)
+! updating velocities
+ velocs_poroelastic(1,iglob) = velocs_poroelastic(1,iglob) + deltatover2*accels_poroelastic(1,iglob)
+ velocs_poroelastic(2,iglob) = velocs_poroelastic(2,iglob) + deltatover2*accels_poroelastic(2,iglob)
+ veloc_elastic(1,iglob) = veloc_elastic(1,iglob) + deltatover2*accel_elastic(1,iglob)
+ veloc_elastic(2,iglob) = veloc_elastic(2,iglob) + deltatover2*accel_elastic(2,iglob)
+! zeros w
+ accelw_poroelastic(1,iglob) = ZERO
+ accelw_poroelastic(2,iglob) = ZERO
+ velocw_poroelastic(1,iglob) = ZERO
+ velocw_poroelastic(2,iglob) = ZERO
+ endif
+
+ enddo
+
+ enddo
+ endif
+
+
!---- compute kinetic and potential energy
if(OUTPUT_ENERGY) &
call compute_energy(displ_elastic,veloc_elastic,displs_poroelastic,velocs_poroelastic, &
Modified: seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt 2009-08-22 15:03:50 UTC (rev 15577)
+++ seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt 2009-08-23 16:28:36 UTC (rev 15578)
@@ -24,15 +24,6 @@
---------------------------
-IMPORTANT KNOWN PROBLEM: the poroelastic solver only works for regular (non
-distorted) meshes, i.e., meshes composed of rectangles aligned with the X and
-Y axes. If the mesh is distorted the solver fails. You can contact
-Christina Morency <cmorency aT princeton DOT edu> for more details.
-
-Dimitri Komatitsch, August 22, 2009.
-
----------------------------
-
- improving compiling with SCOTCH (issue with header file scotchf.h which is Fortran77 legal). Having our own scotchf.h file (without the comments) is not wise.
- comparing the different partitioning methods for METIS and SCOTCH, and finding a good default for SCOTCH.
More information about the CIG-COMMITS
mailing list