[cig-commits] r16412 - in seismo/2D/SPECFEM2D/trunk: . DATA
cmorency at geodynamics.org
cmorency at geodynamics.org
Fri Mar 12 21:38:05 PST 2010
Author: cmorency
Date: 2010-03-12 21:38:05 -0800 (Fri, 12 Mar 2010)
New Revision: 16412
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90
seismo/2D/SPECFEM2D/trunk/checkgrid.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/locate_source_force.F90
seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
- Force source position is now interpolated, similarly to moment tensor source & receiver locations
- For adjoint/kernels calculation the lastframe file is opened at the end of the first time loop leading to a clean reconstructed forward wavefield (Yang Luo)
- Kernels are now written at the local level (i,j,ipsec)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2010-03-13 05:38:05 UTC (rev 16412)
@@ -32,7 +32,7 @@
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# absorbing boundaries parameters
-absorbing_conditions = .true. # absorbing boundary active or not
+absorbing_conditions = .false. # absorbing boundary active or not
absorbbottom = .true.
absorbright = .true.
absorbtop = .false.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2010-03-13 05:38:05 UTC (rev 16412)
@@ -1,20 +1,100 @@
-S0001 AA 75.0000000 -1000.0000000 0.0 0.0
-S0002 AA 75.0000000 -1010.5263158 0.0 0.0
-S0003 AA 75.0000000 -1021.0526316 0.0 0.0
-S0004 AA 75.0000000 -1031.5789474 0.0 0.0
-S0005 AA 75.0000000 -1042.1052632 0.0 0.0
-S0006 AA 75.0000000 -1052.6315789 0.0 0.0
-S0007 AA 75.0000000 -1063.1578947 0.0 0.0
-S0008 AA 75.0000000 -1073.6842105 0.0 0.0
-S0009 AA 75.0000000 -1084.2105263 0.0 0.0
-S0010 AA 75.0000000 -1094.7368421 0.0 0.0
-S0011 AA 75.0000000 -1105.2631579 0.0 0.0
-S0012 AA 75.0000000 -1115.7894737 0.0 0.0
-S0013 AA 75.0000000 -1126.3157895 0.0 0.0
-S0014 AA 75.0000000 -1136.8421053 0.0 0.0
-S0015 AA 75.0000000 -1147.3684211 0.0 0.0
-S0016 AA 75.0000000 -1157.8947368 0.0 0.0
-S0017 AA 75.0000000 -1168.4210526 0.0 0.0
-S0018 AA 75.0000000 -1178.9473684 0.0 0.0
-S0019 AA 75.0000000 -1189.4736842 0.0 0.0
-S0020 AA 75.0000000 -1200.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 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2010-03-13 05:38:05 UTC (rev 16412)
@@ -70,7 +70,7 @@
# GNU gfortran
#F90 = gfortran
-#F90 = mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
+##F90 = mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
#CC = gcc
##FLAGS_NOCHECK = -O3 -march=opteron -m64 -mfpmath=sse,387
#FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O3 -pedantic -pedantic-errors -Wunused -Waliasing -Wampersand -Wline-truncation -Wsurprising -Wunderflow -fno-trapping-math # -mcmodel=medium
@@ -99,7 +99,7 @@
OBJS_SPECFEM2D = $O/checkgrid.o $O/datim.o $O/enforce_acoustic_free_surface.o\
$O/compute_forces_acoustic.o $O/compute_forces_elastic.o\
- $O/compute_forces_solid.o $O/compute_forces_fluid.o\
+ $O/compute_forces_solid.o $O/compute_forces_fluid.o $O/get_poroelastic_velocities.o\
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivation_matrices.o\
$O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o $O/setup_sources_receivers.o\
$O/specfem2D.o $O/write_seismograms.o $O/define_external_model.o $O/createnum_fast.o $O/createnum_slow.o\
@@ -183,6 +183,9 @@
$O/attenuation_model.o: attenuation_model.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/attenuation_model.o attenuation_model.f90
+
+$O/get_poroelastic_velocities.o: get_poroelastic_velocities.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/get_poroelastic_velocities.o get_poroelastic_velocities.f90
### use optimized compilation option for solver only
$O/specfem2D.o: specfem2D.F90 constants.h
Modified: seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2010-03-13 05:38:05 UTC (rev 16412)
@@ -5,7 +5,8 @@
regular and unstructured meshes, generated for example by Cubit
(cubit.sandia.gov). The solver has adjoint capabilities and can
calculate finite-frequency sensitivity kernels for acoustic,
- (an)elastic, and poroelastic media. Finally, the solver can run
+ (an)elastic, and poroelastic media. The package also considers 2D SH
+ and P-SV wave propagation. Finally, the solver can run
both in serial and in parallel. See SPECFEM2D
<http://www.geodynamics.org/cig/software/packages/seismo/specfem2d>
for the source code.
@@ -272,6 +273,12 @@
Note: for the poroelastic case, mu_s is irrelevant.
For details on the poroelastic theory see Morency and Tromp, GJI 2008.
+get_poroelastic_velocities.f90 allows to compute cpI, cpII, and cs function of
+the source dominant frequency. Notice that for this calculation we use permxx
+and the dominant frequency of the first source , f0(1). Caution if you use
+several sources with different frequencies and if you consider anistropic
+permeability.
+
--------------------------------------------------
HOW TO OBTAIN FINITE SENSITIVITY KERNELS
--------------------------------------------------
@@ -313,7 +320,8 @@
Output_files (for example for the elastic case)
snapshot_rho_kappa_mu*****
snapshot_rhop_alpha_beta*****
-which are the primary moduli kernels and the phase velocities kernels respectively.
+which are the primary moduli kernels and the phase velocities kernels respectively, in ascii format
+and at the local level, that is as "kernels(i,j,ispec)".
Note1: At the moment, adjoint simulations do not support anisotropy, attenuation, and viscous damping.
Note2: You will need S****.AA.BHX.adj, S****.AA.BHY.adj and S****.AA.BHZ.adj
Modified: seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/adj_seismogram.f90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -6,10 +6,10 @@
implicit none
!
!!!! user edit
- integer, parameter :: NSTEP = 26000
+ integer, parameter :: NSTEP = 3000
integer, parameter :: nrec = 1
- double precision, parameter :: t0 = 8
- double precision, parameter :: deltat = 0.25d-2
+ double precision, parameter :: t0 = 12
+ double precision, parameter :: deltat = 6d-2
double precision, parameter :: EPS = 1.d-40
!!!!
integer :: itime,icomp,istart,iend,nlen,irec,NDIM,NDIMr,adj_comp
@@ -26,17 +26,17 @@
!!!! user edit
! which calculation: P-SV (use (1)) or SH (membrane) (use (2)) waves
-! NDIMr=2 !(1)
- NDIMr=1 !(2)
+ NDIMr=2 !(1)
+! NDIMr=1 !(2)
! list of stations
station_name(1) = 'S0001'
- tstart(1) = 39.7d0 + t0
- tend(1) = 41d0 + t0
+ tstart(1) = 100d0 + t0
+ tend(1) = 120d0 + t0
! which calculation: P-SV (use (1)) or SH (membrane) (use (2)) waves
-! compr = (/"BHX","BHZ"/) !(1)
- compr = (/"BHY","dummy"/) !(2)
+ compr = (/"BHX","BHZ"/) !(1)
+! compr = (/"BHY","dummy"/) !(2)
! chose the component for the adjoint source (adj_comp = 1: X, 2:Y, 3:Z)
- adj_comp = 2
+ adj_comp = 1
!!!!
do irec =1,nrec
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -42,10 +42,11 @@
!
!========================================================================
- subroutine checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,ibool,kmato,coord,npoin, &
+ subroutine checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,permeability,ibool,kmato,coord,npoin, &
vpImin,vpImax,vpIImin,vpIImax,assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat, &
f0,t0,initialfield,time_function_type,coorg,xinterp,zinterp,shapeint,knods,simulation_title, &
- npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic,myrank,nproc,NSOURCE,poroelastic)
+ npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic,myrank,nproc,NSOURCE,poroelastic, &
+ freq0,Q0,TURN_VISCATTENUATION_ON)
! check the mesh, stability and number of points per wavelength
@@ -79,6 +80,7 @@
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: poroelastcoef
double precision, dimension(numat) :: porosity,tortuosity
+ double precision, dimension(3,numat) :: permeability
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
double precision coord(NDIM,npoin)
@@ -86,15 +88,16 @@
double precision vpImin,vpImax,vsmin,vsmax,densmin,densmax,vpImax_local,vpImin_local,vsmin_local
double precision vpIImin,vpIImax,vpIImax_local,vpIImin_local
double precision kappa_s,kappa_f,kappa_fr,mu_s,mu_fr,denst_s,denst_f,denst,phi,tort,cpIloc,cpIIloc,csloc
- double precision afactor,bfactor,cfactor,D_biot,H_biot,C_biot,M_biot,cpIsquare,cpIIsquare,cssquare
- double precision f0min,f0max
+ double precision D_biot,H_biot,C_biot,M_biot,cpIsquare,cpIIsquare,cssquare
+ double precision f0min,f0max,freq0,Q0,w_c,eta_f,perm
double precision lambdaplus2mu,mu
double precision distance_min,distance_max,distance_min_local,distance_max_local
double precision courant_stability_number_max,lambdaPImin,lambdaPImax,lambdaPIImin,lambdaPIImax, &
lambdaSmin,lambdaSmax
double precision deltat,distance_1,distance_2,distance_3,distance_4
double precision, dimension(NSOURCE) :: f0,t0
- logical assign_external_model,initialfield,any_elastic,any_poroelastic,all_anisotropic
+ logical assign_external_model,initialfield,any_elastic,any_poroelastic,all_anisotropic, &
+ TURN_VISCATTENUATION_ON
! for the stability condition
! maximum polynomial degree for which we can compute the stability condition
@@ -1417,6 +1420,7 @@
if(poroelastic(ispec)) then
phi = porosity(material)
tort = tortuosity(material)
+ perm = permeability(1,material)
!solid properties
mu_s = poroelastcoef(2,1,material)
kappa_s = poroelastcoef(3,1,material) - FOUR_THIRDS*mu_s
@@ -1424,23 +1428,19 @@
!fluid properties
kappa_f = poroelastcoef(1,2,material)
denst_f = density(2,material)
+ eta_f = poroelastcoef(2,2,material)
!frame properties
mu_fr = poroelastcoef(2,3,material)
kappa_fr = poroelastcoef(3,3,material) - FOUR_THIRDS*mu_fr
- denst = (1.d0 - phi)*denst_s + phi*denst_f
!Biot coefficients for the input phi
D_biot = kappa_s*(1.d0 + phi*(kappa_s/kappa_f - 1.d0))
H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
-! Approximated velocities (no viscous dissipation)
- afactor = denst - phi/tort*denst_f
- bfactor = H_biot + phi*denst/(tort*denst_f)*M_biot - 2.d0*phi/tort*C_biot
- cfactor = phi/(tort*denst_f)*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cssquare = mu_fr/afactor
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mu_fr,phi, &
+ tort,denst_s,denst_f,eta_f,perm,f0(1),freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
cpIloc = sqrt(cpIsquare)
cpIIloc = sqrt(cpIIsquare)
csloc = sqrt(cssquare)
@@ -1881,6 +1881,7 @@
if(poroelastic(ispec)) then
phi=porosity(material)
tort=tortuosity(material)
+ perm=permeability(1,material)
!solid properties
mu_s = poroelastcoef(2,1,material)
kappa_s = poroelastcoef(3,1,material) - FOUR_THIRDS*mu_s
@@ -1888,21 +1889,19 @@
!fluid properties
kappa_f = poroelastcoef(1,2,material)
denst_f = density(2,material)
+ eta_f = poroelastcoef(2,2,material)
!frame properties
mu_fr = poroelastcoef(2,3,material)
kappa_fr = poroelastcoef(3,3,material) - FOUR_THIRDS*mu_fr
- denst = (1.d0 - phi)*denst_s + phi*denst_f
!Biot coefficients for the input phi
D_biot = kappa_s*(1.d0 + phi*(kappa_s/kappa_f - 1.d0))
H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
-! Approximated velocities (no viscous dissipation)
- afactor = denst - phi/tort*denst_f
- bfactor = H_biot + phi*denst/(tort*denst_f)*M_biot - 2.d0*phi/tort*C_biot
- cfactor = phi/(tort*denst_f)*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mu_fr,phi, &
+ tort,denst_s,denst_f,eta_f,perm,f0(1),freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
cpIloc = sqrt(cpIsquare)
else
lambdaplus2mu = poroelastcoef(3,1,material)
@@ -2233,6 +2232,7 @@
if(poroelastic(ispec)) then
phi = porosity(material)
tort = tortuosity(material)
+ perm = permeability(1,material)
!solid properties
mu_s = poroelastcoef(2,1,material)
kappa_s = poroelastcoef(3,1,material) - FOUR_THIRDS*mu_s
@@ -2240,23 +2240,19 @@
!fluid properties
kappa_f = poroelastcoef(1,2,material)
denst_f = density(2,material)
+ eta_f = poroelastcoef(2,2,material)
!frame properties
mu_fr = poroelastcoef(2,3,material)
kappa_fr = poroelastcoef(3,3,material) - FOUR_THIRDS*mu_fr
- denst = (1.d0 - phi)*denst_s + phi*denst_f
!Biot coefficients for the input phi
D_biot = kappa_s*(1.d0 + phi*(kappa_s/kappa_f - 1.d0))
H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
-! Approximated velocities (no viscous dissipation)
- afactor = denst - phi/tort*denst_f
- bfactor = H_biot + phi*denst/(tort*denst_f)*M_biot - 2.d0*phi/tort*C_biot
- cfactor = phi/(tort*denst_f)*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cssquare = mu_fr/afactor
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mu_fr,phi, &
+ tort,denst_s,denst_f,eta_f,perm,f0(1),freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
cpIloc = sqrt(cpIsquare)
csloc = sqrt(cssquare)
else
@@ -2651,6 +2647,7 @@
if(poroelastic(ispec)) then
phi = porosity(material)
tort = tortuosity(material)
+ perm = permeability(1,material)
!solid properties
mu_s = poroelastcoef(2,1,material)
kappa_s = poroelastcoef(3,1,material) - FOUR_THIRDS*mu_s
@@ -2658,20 +2655,19 @@
!fluid properties
kappa_f = poroelastcoef(1,2,material)
denst_f = density(2,material)
+ eta_f = poroelastcoef(2,2,material)
!frame properties
mu_fr = poroelastcoef(2,3,material)
kappa_fr = poroelastcoef(3,3,material) - FOUR_THIRDS*mu_fr
- denst = (1.d0 - phi)*denst_s + phi*denst_f
!Biot coefficients for the input phi
D_biot = kappa_s*(1.d0 + phi*(kappa_s/kappa_f - 1.d0))
H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
-! Approximated velocities (no viscous dissipation)
- afactor = denst - phi/tort*denst_f
- bfactor = H_biot + phi*denst/(tort*denst_f)*M_biot - 2.d0*phi/tort*C_biot
- cfactor = phi/(tort*denst_f)*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
+
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mu_fr,phi, &
+ tort,denst_s,denst_f,eta_f,perm,f0(1),freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
cpIloc = sqrt(cpIsquare)
else
lambdaplus2mu = poroelastcoef(3,1,material)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_fluid.f90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -60,7 +60,7 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
C_k,M_k,NSOURCE,nrec,isolver,save_forward,&
b_absorb_poro_w_left,b_absorb_poro_w_right,b_absorb_poro_w_bottom,b_absorb_poro_w_top,&
- nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax)
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,f0,freq0,Q0)
! compute forces for the fluid poroelastic part
@@ -134,7 +134,10 @@
real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: wzgll
+!
+ double precision :: f0,freq0,Q0,w_c
+
!---
!--- local variables
!---
@@ -171,7 +174,7 @@
real(kind=CUSTOM_REAL) :: etal_f,kappal_f,rhol_f
real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr,phil,tortl,viscodampx,viscodampz
real(kind=CUSTOM_REAL) :: permlxx,permlxz,permlzz,invpermlxx,invpermlxz,invpermlzz,detk
- real(kind=CUSTOM_REAL) :: afactor,bfactor,cfactor,D_biot,H_biot,C_biot,M_biot,rhol_bar
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,rhol_bar
real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,cpIl,cpIIl,csl
@@ -528,6 +531,7 @@
! get poroelastic parameters of current spectral element
phil = porosity(kmato(ispec))
tortl = tortuosity(kmato(ispec))
+ permlxx = permeability(1,kmato(ispec))
!solid properties
mul_s = poroelastcoef(2,1,kmato(ispec))
kappal_s = poroelastcoef(3,1,kmato(ispec)) - 4._CUSTOM_REAL*mul_s/3._CUSTOM_REAL
@@ -535,6 +539,7 @@
!fluid properties
kappal_f = poroelastcoef(1,2,kmato(ispec))
rhol_f = density(2,kmato(ispec))
+ etal_f = poroelastcoef(2,2,kmato(ispec))
!frame properties
mul_fr = poroelastcoef(2,3,kmato(ispec))
kappal_fr = poroelastcoef(3,3,kmato(ispec)) - 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
@@ -544,14 +549,10 @@
H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + kappal_fr + 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
-! Approximated velocities (no viscous dissipation)
- afactor = rhol_bar - phil/tortl*rhol_f
- bfactor = H_biot + phil*rhol_bar/(tortl*rhol_f)*M_biot - TWO*phil/tortl*C_biot
- cfactor = phil/(tortl*rhol_f)*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(TWO*afactor)
- cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(TWO*afactor)
- cssquare = mul_fr/afactor
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mul_fr,phil, &
+ tortl,rhol_s,rhol_f,etal_f,permlxx,f0,freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
cpIl = sqrt(cpIsquare)
cpIIl = sqrt(cpIIsquare)
csl = sqrt(cssquare)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_solid.f90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -60,7 +60,7 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
mufr_k,B_k,NSOURCE,nrec,isolver,save_forward,&
b_absorb_poro_s_left,b_absorb_poro_s_right,b_absorb_poro_s_bottom,b_absorb_poro_s_top,&
- nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax)
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,f0,freq0,Q0)
! compute forces for the solid poroelastic part
@@ -134,6 +134,8 @@
real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: wzgll
+!
+ double precision :: f0,freq0,Q0,w_c
!---
!--- local variables
@@ -174,7 +176,7 @@
real(kind=CUSTOM_REAL) :: etal_f,kappal_f,rhol_f
real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr,phil,tortl,viscodampx,viscodampz
real(kind=CUSTOM_REAL) :: permlxx,permlxz,permlzz,invpermlxx,invpermlxz,invpermlzz,detk
- real(kind=CUSTOM_REAL) :: afactor,bfactor,cfactor,D_biot,H_biot,C_biot,M_biot,rhol_bar
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,rhol_bar
real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,cpIl,cpIIl,csl
@@ -547,6 +549,7 @@
! get poroelastic parameters of current spectral element
phil = porosity(kmato(ispec))
tortl = tortuosity(kmato(ispec))
+ permlxx = permeability(1,kmato(ispec))
!solid properties
mul_s = poroelastcoef(2,1,kmato(ispec))
kappal_s = poroelastcoef(3,1,kmato(ispec)) - 4._CUSTOM_REAL*mul_s/3._CUSTOM_REAL
@@ -554,6 +557,7 @@
!fluid properties
kappal_f = poroelastcoef(1,2,kmato(ispec))
rhol_f = density(2,kmato(ispec))
+ etal_f = poroelastcoef(2,2,kmato(ispec))
!frame properties
mul_fr = poroelastcoef(2,3,kmato(ispec))
kappal_fr = poroelastcoef(3,3,kmato(ispec)) - 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
@@ -563,14 +567,10 @@
H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + kappal_fr + 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
-! Approximated velocities (no viscous dissipation)
- afactor = rhol_bar - phil/tortl*rhol_f
- bfactor = H_biot + phil*rhol_bar/(tortl*rhol_f)*M_biot - TWO*phil/tortl*C_biot
- cfactor = phil/(tortl*rhol_f)*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(TWO*afactor)
- cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(TWO*afactor)
- cssquare = mul_fr/afactor
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mul_fr,phil, &
+ tortl,rhol_s,rhol_f,etal_f,permlxx,f0,freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
cpIl = sqrt(cpIsquare)
cpIIl = sqrt(cpIIsquare)
csl = sqrt(cssquare)
Modified: seismo/2D/SPECFEM2D/trunk/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gmat01.f90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/gmat01.f90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -42,8 +42,8 @@
!
!========================================================================
-subroutine gmat01(density_array,porosity_array,tortuosity_array,aniso_array,permeability,poroelastcoef,&
- numat,myrank,ipass,Qp_array,Qs_array)
+ subroutine gmat01(density_array,porosity_array,tortuosity_array,aniso_array,permeability,poroelastcoef,&
+ numat,myrank,ipass,Qp_array,Qs_array,freq0,Q0,f0,TURN_VISCATTENUATION_ON)
! read properties of a 2D isotropic or anisotropic linear elastic element
@@ -67,7 +67,9 @@
double precision val1,val2,val3,val4,val5,val6
double precision val7,val8,val9,val10,val11,val12,val0
double precision :: c11,c13,c15,c33,c35,c55
- double precision afactor,bfactor,cfactor,D_biot,H_biot,C_biot,M_biot,density_bar
+ double precision D_biot,H_biot,C_biot,M_biot
+ double precision f0,Q0,freq0,w_c
+ logical TURN_VISCATTENUATION_ON
!
!---- loop over the different material sets
@@ -190,15 +192,10 @@
H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
- ! Approximated velocities (no viscous dissipation)
- density_bar = (1.d0 - phi)*density(1) + phi*density(2)
- afactor = density_bar - phi/tortuosity*density(2)
- bfactor = H_biot + phi*density_bar/(tortuosity*density(2))*M_biot - 2.d0*phi/tortuosity*C_biot
- cfactor = phi/(tortuosity*density(2))*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cssquare = val11/afactor
+ call get_poroelastic_velocities(cpIsquare,cpIIsquare,cssquare,H_biot,C_biot,M_biot,mu_fr,phi, &
+ tortuosity,density(1),density(2),eta_f,val4,f0,freq0,Q0,w_c,TURN_VISCATTENUATION_ON)
+
porosity_array(n) = val2
tortuosity_array(n) = val3
permeability(1,n) = val4
@@ -296,7 +293,7 @@
write(iout,700) density(2),kappa_f,eta_f
write(iout,800) lambda_fr,mu_fr,kappa_fr,porosity_array(n),tortuosity_array(n),&
permeability(1,n),permeability(2,n),permeability(3,n),Qs
- write(iout,900) D_biot,H_biot,C_biot,M_biot
+ write(iout,900) D_biot,H_biot,C_biot,M_biot,w_c
endif
endif
@@ -381,13 +378,14 @@
'Permeability zz component. . . . . . . . . . =',1pe15.8,/5x,&
'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
-900 format(//5x,'-------------------------------',/5x, &
- '-- Biot coefficients --',/5x, &
- '-------------------------------',/5x, &
- 'D. . . . . . . . =',1pe15.8,/5x, &
- 'H. . . . . . . . =',1pe15.8,/5x, &
- 'C. . . . . . . . =',1pe15.8,/5x, &
- 'M. . . . . . . . =',1pe15.8)
+900 format(//5x,'-------------------------------',/5x, &
+ '-- Biot coefficients --',/5x, &
+ '-------------------------------',/5x, &
+ 'D. . . . . . . . =',1pe15.8,/5x, &
+ 'H. . . . . . . . =',1pe15.8,/5x, &
+ 'C. . . . . . . . =',1pe15.8,/5x, &
+ 'M. . . . . . . . =',1pe15.8,/5x, &
+ 'characteristic freq =',1pe15.8)
end subroutine gmat01
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -42,12 +42,13 @@
!
!========================================================================
- subroutine locate_source_force(coord,ibool,npoin,nspec,x_source,z_source,ix_source,iz_source, &
- ispec_source,iglob_source,is_proc_source,nb_proc_source,ipass)
+!----
+!---- locate_source_force finds the correct position of the point force source
+!----
-!
-!----- calculer la position reelle de la source
-!
+ subroutine locate_source_force(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
+ ispec_selected_source,is_proc_source,nb_proc_source,nproc,myrank, &
+ xi_source,gamma_source,coorg,knods,ngnod,npgeo,ipass)
implicit none
@@ -56,101 +57,187 @@
include "mpif.h"
#endif
- integer npoin,nspec,ipass
- integer ibool(NGLLX,NGLLZ,nspec)
+ integer nspec,npoin,ngnod,npgeo,ipass
- double precision x_source,z_source
+ integer knods(ngnod,nspec)
+ double precision coorg(NDIM,npgeo)
+
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+
+! array containing coordinates of the points
double precision coord(NDIM,npoin)
- integer, intent(inout) :: is_proc_source, nb_proc_source
+ integer i,j,ispec,iglob,iter_loop,ix_initial_guess,iz_initial_guess
- integer ip,ix,iz,numelem,ilowx,ilowz,ihighx,ihighz,ix_source,iz_source,ispec_source,iglob_source
+ double precision x_source,z_source,dist
+ double precision xi,gamma,dx,dz,dxi,dgamma
- double precision distminmax,distmin,xp,zp,dist, dist_glob
+! Gauss-Lobatto-Legendre points of integration
+ double precision xigll(NGLLX)
+ double precision zigll(NGLLZ)
+ double precision x,z,xix,xiz,gammax,gammaz,jacobian
+ double precision distmin,final_distance,dist_glob
+
+! source information
+ integer ispec_selected_source,is_proc_source,nb_proc_source
+ integer, intent(in) :: nproc, myrank
+ double precision xi_source,gamma_source
+
+#ifdef USE_MPI
+ integer, dimension(1:nproc) :: allgather_is_proc_source
+ integer, dimension(1) :: locate_is_proc_source
integer :: ierror
+#endif
- ierror = 0
- is_proc_source = 0
- distminmax = -HUGEVAL
- distmin = +HUGEVAL
+! **************
+ if ((myrank == 0 .or. nproc == 1) .and. ipass == 1) then
+ write(IOUT,*)
+ write(IOUT,*) '*******************************'
+ write(IOUT,*) ' locating force source'
+ write(IOUT,*) '*******************************'
+ write(IOUT,*)
+ endif
- ilowx = 1
- ilowz = 1
- ihighx = NGLLX
- ihighz = NGLLZ
+! set distance to huge initial value
+ distmin = HUGEVAL
-! look for the closest grid point
- do numelem = 1,nspec
+ is_proc_source = 0
- do ix = ilowx,ihighx
- do iz = ilowz,ihighz
+ do ispec = 1,nspec
-! global point number
- ip = ibool(ix,iz,numelem)
+! loop only on points inside the element
+! exclude edges to ensure this point is not shared with other elements
+ do j = 2,NGLLZ-1
+ do i = 2,NGLLX-1
-! coordinates of this grid point
- xp = coord(1,ip)
- zp = coord(2,ip)
+ iglob = ibool(i,j,ispec)
+ dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
- dist = sqrt((xp-x_source)**2 + (zp-z_source)**2)
-
-! keep the point for which distance is minimum
- if(dist < distmin) then
+! keep this point if it is closer to the source
+ if(dist < distmin) then
distmin = dist
- iglob_source = ip
- ix_source = ix
- iz_source = iz
- ispec_source = numelem
- endif
+ ispec_selected_source = ispec
+ ix_initial_guess = i
+ iz_initial_guess = j
+ endif
- enddo
- enddo
+ enddo
+ enddo
- enddo
+! end of loop on all the spectral elements
+ enddo
- distminmax = max(distmin,distminmax)
-
#ifdef USE_MPI
-! global minimum distance computed over all processes
- call MPI_ALLREDUCE (distminmax, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
+ ! global minimum distance computed over all processes
+ call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
#else
- dist_glob = distminmax
+ dist_glob = distmin
#endif
! check if this process contains the source
- if (dist_glob == distminmax) is_proc_source = 1
+ if ( dist_glob == distmin ) is_proc_source = 1
#ifdef USE_MPI
-! determining the number of processes that contain the source (useful when the source is located on an interface)
+ ! determining the number of processes that contain the source (useful when the source is located on an interface)
call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
#else
nb_proc_source = is_proc_source
+
#endif
- if (nb_proc_source < 1) call exit_MPI('error locating force source')
+#ifdef USE_MPI
+ ! when several processes contain the source, we elect one of them (minimum rank).
+ if ( nb_proc_source > 1 ) then
+
+ call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
+ locate_is_proc_source = maxloc(allgather_is_proc_source) - 1
+
+ if ( myrank /= locate_is_proc_source(1) ) then
+ is_proc_source = 0
+ endif
+ nb_proc_source = 1
+
+ endif
+
+#endif
+
+! ****************************************
+! find the best (xi,gamma) for each source
+! ****************************************
+
+! use initial guess in xi and gamma
+ xi = xigll(ix_initial_guess)
+ gamma = zigll(iz_initial_guess)
+
+! iterate to solve the non linear system
+ do iter_loop = 1,NUM_ITER
+
+! recompute jacobian for the new point
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo)
+
+! compute distance to target location
+ dx = - (x - x_source)
+ dz = - (z - z_source)
+
+! compute increments
+ dxi = xix*dx + xiz*dz
+ dgamma = gammax*dx + gammaz*dz
+
+! update values
+ xi = xi + dxi
+ gamma = gamma + dgamma
+
+! impose that we stay in that element
+! (useful if user gives a source outside the mesh for instance)
+! we can go slightly outside the [1,1] segment since with finite elements
+! the polynomial solution is defined everywhere
+! this can be useful for convergence of itertive scheme with distorted elements
+ if (xi > 1.10d0) xi = 1.10d0
+ if (xi < -1.10d0) xi = -1.10d0
+ if (gamma > 1.10d0) gamma = 1.10d0
+ if (gamma < -1.10d0) gamma = -1.10d0
+
+! end of non linear iterations
+ enddo
+
+! compute final coordinates of point found
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo)
+
+! store xi,gamma of point found
+ xi_source = xi
+ gamma_source = gamma
+
+! compute final distance between asked and found
+ final_distance = sqrt((x_source-x)**2 + (z_source-z)**2)
+
if (is_proc_source == 1 .and. ipass == 1) then
- write(IOUT,200)
- write(IOUT,"(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3,1x,i5.5)") x_source,z_source, &
- coord(1,iglob_source),coord(2,iglob_source),distmin,nb_proc_source
write(IOUT,*)
+ write(IOUT,*) 'Force source:'
+
+ if(final_distance == HUGEVAL) call exit_MPI('error locating force source')
+
+ write(IOUT,*) ' original x: ',sngl(x_source)
+ write(IOUT,*) ' original z: ',sngl(z_source)
+ write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
+ write(IOUT,*) ' in element ',ispec_selected_source
+ write(IOUT,*) ' at xi,gamma coordinates = ',xi_source,gamma_source
write(IOUT,*)
- write(IOUT,"('Maximum distance between asked and real =',f12.3)") distminmax
+
+ write(IOUT,*)
+ write(IOUT,*) 'end of force source detection'
+ write(IOUT,*)
endif
#ifdef USE_MPI
call MPI_BARRIER(MPI_COMM_WORLD,ierror)
#endif
- 200 format(//1x,48('=')/,' = S o u r c e s ', &
- 'r e a l p o s i t i o n s ='/1x,48('=')// &
- ' Source x-asked z-asked x-obtain z-obtain dist nb_proc_source'/)
-
end subroutine locate_source_force
Modified: seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -43,11 +43,11 @@
!
!========================================================================
-subroutine setup_sources_receivers(NSOURCE,assign_external_model,initialfield,numat,source_type,&
- coord,ibool,kmato,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
- x_source,z_source,ix_source,iz_source,ispec_selected_source,ispec_selected_rec,iglob_source, &
- is_proc_source,nb_proc_source,rho_at_source_location,ipass,&
- sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,density,rhoext,&
+subroutine setup_sources_receivers(NSOURCE,initialfield,source_type,&
+ coord,ibool,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
+ x_source,z_source,ispec_selected_source,ispec_selected_rec, &
+ is_proc_source,nb_proc_source,ipass,&
+ sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,&
nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod, &
nrec,nrecloc,recloc,which_proc_receiver,st_xval,st_zval, &
xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver)
@@ -56,8 +56,8 @@
include "constants.h"
- logical :: assign_external_model,initialfield
- integer :: NSOURCE, numat
+ logical :: initialfield
+ integer :: NSOURCE
integer :: npgeo,ngnod,myrank,ipass,nproc
integer :: npoin,nspec,nelem_acoustic_surface
@@ -77,63 +77,64 @@
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
! for sources
- double precision :: rho_at_source_location
integer, dimension(NSOURCE) :: source_type
- integer, dimension(NSOURCE) :: ispec_selected_source,iglob_source,ix_source,iz_source,is_proc_source,nb_proc_source
+ integer, dimension(NSOURCE) :: ispec_selected_source,is_proc_source,nb_proc_source
real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
double precision, dimension(NSOURCE) :: x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz
logical, dimension(nspec) :: elastic,poroelastic
- integer, dimension(nspec) :: kmato
integer, dimension(ngnod,nspec) :: knods
integer, dimension(5,nelem_acoustic_surface) :: acoustic_surface
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
double precision, dimension(NDIM,npgeo) :: coorg
- double precision, dimension(2,numat) :: density
- double precision, dimension(NGLLX,NGLLZ,nspec) :: rhoext
double precision, dimension(NDIM,npoin) :: coord
+ integer :: ixmin, ixmax, izmin, izmax
+
! Local variables
- integer i_source,ispec,ispec_acoustic_surface,i,j,iglob
+ integer i_source,ispec,ispec_acoustic_surface
do i_source=1,NSOURCE
if(source_type(i_source) == 1) then
! collocated force source
- call locate_source_force(coord,ibool,npoin,nspec,x_source(i_source),z_source(i_source), &
- ix_source(i_source),iz_source(i_source),ispec_selected_source(i_source),iglob_source(i_source), &
- is_proc_source(i_source),nb_proc_source(i_source),ipass)
+ call locate_source_force(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
+ ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),&
+ nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass)
- ! get density at the source in order to implement collocated force with the right
- ! amplitude later
- if(is_proc_source(i_source) == 1) then
- rho_at_source_location = density(1,kmato(ispec_selected_source(i_source)))
- ! external velocity model
- if(assign_external_model) rho_at_source_location = &
- rhoext(ix_source(i_source),iz_source(i_source),ispec_selected_source(i_source))
- endif
+! check that acoustic source is not exactly on the free surface because pressure is zero there
+ if(is_proc_source(i_source) == 1) then
+ do ispec_acoustic_surface = 1,nelem_acoustic_surface
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+ ixmin = acoustic_surface(2,ispec_acoustic_surface)
+ ixmax = acoustic_surface(3,ispec_acoustic_surface)
+ izmin = acoustic_surface(4,ispec_acoustic_surface)
+ izmax = acoustic_surface(5,ispec_acoustic_surface)
+ if( .not. elastic(ispec) .and. .not. poroelastic(ispec) .and. ispec == ispec_selected_source(i_source) ) then
+ if ( (izmin==1 .and. izmax==1 .and. ixmin==1 .and. ixmax==NGLLX .and. &
+ gamma_source(i_source) < -0.99d0) .or.&
+ (izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==NGLLX .and. &
+ gamma_source(i_source) > 0.99d0) .or.&
+ (izmin==1 .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==1 .and. &
+ xi_source(i_source) < -0.99d0) .or.&
+ (izmin==1 .and. izmax==NGLLZ .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
+ xi_source(i_source) > 0.99d0) .or.&
+ (izmin==1 .and. izmax==1 .and. ixmin==1 .and. ixmax==1 .and. &
+ gamma_source(i_source) < -0.99d0 .and. xi_source(i_source) < -0.99d0) .or.&
+ (izmin==1 .and. izmax==1 .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
+ gamma_source(i_source) < -0.99d0 .and. xi_source(i_source) > 0.99d0) .or.&
+ (izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==1 .and. &
+ gamma_source(i_source) > 0.99d0 .and. xi_source(i_source) < -0.99d0) .or.&
+ (izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
+ gamma_source(i_source) > 0.99d0 .and. xi_source(i_source) > 0.99d0) ) then
+ call exit_MPI('an acoustic source cannot be located exactly on the free surface because pressure is zero there')
+ endif
+ endif
+ enddo
+ endif
- ! check that acoustic source is not exactly on the free surface because pressure is zero there
- if(is_proc_source(i_source) == 1) then
- do ispec_acoustic_surface = 1,nelem_acoustic_surface
- ispec = acoustic_surface(1,ispec_acoustic_surface)
- if( .not. elastic(ispec) .and. .not. poroelastic(ispec) .and. ispec == ispec_selected_source(i_source) ) then
- do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
- do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
- iglob = ibool(i,j,ispec)
- if ( iglob_source(i_source) == iglob ) then
-
- call exit_MPI('an acoustic source cannot be located exactly on the free surface because pressure is zero there')
-
- endif
- enddo
- enddo
- endif
- enddo
- endif
-
else if(source_type(i_source) == 2) then
! moment-tensor source
call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2010-03-13 02:23:59 UTC (rev 16411)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2010-03-13 05:38:05 UTC (rev 16412)
@@ -316,7 +316,6 @@
double precision, dimension(:), allocatable :: vp_display
double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
- double precision :: rho_at_source_location
double precision, dimension(:,:,:), allocatable :: Qp_attenuationext,Qs_attenuationext
double precision, dimension(:,:,:), allocatable :: c11ext,c13ext,c15ext,c33ext,c35ext,c55ext
@@ -330,7 +329,7 @@
integer, dimension(:), allocatable :: kmato,numabs, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top,jbegin_left,jend_left,jbegin_right,jend_right
- integer, dimension(:), allocatable :: ispec_selected_source,iglob_source,ix_source,iz_source,&
+ integer, dimension(:), allocatable :: ispec_selected_source,iglob_source,&
is_proc_source,nb_proc_source
double precision displnorm_all,displnorm_all_glob
double precision, dimension(:), allocatable :: aval
@@ -438,14 +437,14 @@
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accel_elastic,b_veloc_elastic,b_displ_elastic
real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_potential_dot_dot_acoustic,b_potential_dot_acoustic,b_potential_acoustic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_ac,b_displ_ac,b_accel_ac
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_kl, mu_kl, kappa_kl
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rho_kl, mu_kl, kappa_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhol_global, mul_global, kappal_global
real(kind=CUSTOM_REAL), dimension(:), allocatable :: mu_k, kappa_k,rho_k
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhop_kl, beta_kl, alpha_kl
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_ac_kl, kappa_ac_kl
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhop_kl, beta_kl, alpha_kl
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rho_ac_kl, kappa_ac_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhol_ac_global, kappal_ac_global
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhop_ac_kl, alpha_ac_kl
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhot_kl, rhof_kl, sm_kl, eta_kl, mufr_kl, B_kl, &
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhop_ac_kl, alpha_ac_kl
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhot_kl, rhof_kl, sm_kl, eta_kl, mufr_kl, B_kl, &
C_kl, M_kl, rhob_kl, rhofb_kl, phi_kl, Bb_kl, Cb_kl, Mb_kl, mufrb_kl, &
rhobb_kl, rhofbb_kl, phib_kl, cpI_kl, cpII_kl, cs_kl, ratio_kl
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhot_k, rhof_k, sm_k, eta_k, mufr_k, B_k, &
@@ -504,6 +503,10 @@
double precision, dimension(:), allocatable :: hxir,hgammar,hpxir,hpgammar
double precision, dimension(:,:), allocatable :: hxir_store,hgammar_store
+! Lagrange interpolators at sources
+ double precision, dimension(:), allocatable :: hxis,hgammas,hpxis,hpgammas
+ double precision, dimension(:,:), allocatable :: hxis_store,hgammas_store
+
! for Lagrange interpolants
double precision, external :: hgll
@@ -639,9 +642,9 @@
!!!!!!!!!!
double precision, dimension(:,:,:),allocatable:: rho_local,vp_local,vs_local
!!!! hessian
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhorho_el_hessian_final1, rhorho_el_hessian_temp1
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhorho_el_hessian_final2, rhorho_el_hessian_temp2
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhorho_ac_hessian_final1, rhorho_ac_hessian_final2
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhorho_el_hessian_final1, rhorho_el_hessian_final2
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhorho_el_hessian_temp1, rhorho_el_hessian_temp2
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhorho_ac_hessian_final1, rhorho_ac_hessian_final2
real(kind=CUSTOM_REAL), dimension(:), allocatable :: weight_line_x, weight_line_z, weight_surface,weight_jacobian
integer, dimension(:), allocatable :: weight_gll
real(kind=CUSTOM_REAL) :: zmin_yang, zmax_yang, xmin_yang, xmax_yang
@@ -820,8 +823,6 @@
allocate( ispec_selected_source(NSOURCE) )
allocate( iglob_source(NSOURCE) )
allocate( source_courbe_eros(NSOURCE) )
- allocate( ix_source(NSOURCE) )
- allocate( iz_source(NSOURCE) )
allocate( xi_source(NSOURCE) )
allocate( gamma_source(NSOURCE) )
allocate( is_proc_source(NSOURCE) )
@@ -998,7 +999,7 @@
!---- read the material properties
!
call gmat01(density,porosity,tortuosity,anisotropy,permeability,poroelastcoef,numat,&
- myrank,ipass,Qp_attenuation,Qs_attenuation)
+ myrank,ipass,Qp_attenuation,Qs_attenuation,freq0,Q0,f0(1),TURN_VISCATTENUATION_ON)
!
!---- read spectral macrobloc data
!
@@ -1598,14 +1599,22 @@
! allocate 1-D Lagrange interpolators and derivatives
allocate(hxir(NGLLX))
+ allocate(hxis(NGLLX))
allocate(hpxir(NGLLX))
+ allocate(hpxis(NGLLX))
allocate(hgammar(NGLLZ))
+ allocate(hgammas(NGLLZ))
allocate(hpgammar(NGLLZ))
+ allocate(hpgammas(NGLLZ))
! allocate Lagrange interpolators for receivers
allocate(hxir_store(nrec,NGLLX))
allocate(hgammar_store(nrec,NGLLZ))
+! allocate Lagrange interpolators for sources
+ allocate(hxis_store(NSOURCE,NGLLX))
+ allocate(hgammas_store(NSOURCE,NGLLZ))
+
! allocate other global arrays
allocate(coord(NDIM,npoin))
@@ -1784,11 +1793,11 @@
!---- define actual location of source and receivers
-call setup_sources_receivers(NSOURCE,assign_external_model,initialfield,numat,source_type,&
- coord,ibool,kmato,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
- x_source,z_source,ix_source,iz_source,ispec_selected_source,ispec_selected_rec,iglob_source, &
- is_proc_source,nb_proc_source,rho_at_source_location,ipass,&
- sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,density,rhoext,&
+call setup_sources_receivers(NSOURCE,initialfield,source_type,&
+ coord,ibool,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
+ x_source,z_source,ispec_selected_source,ispec_selected_rec, &
+ is_proc_source,nb_proc_source,ipass,&
+ sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,&
nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod, &
nrec,nrecloc,recloc,which_proc_receiver,st_xval,st_zval, &
xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver)
@@ -2092,6 +2101,14 @@
hgammar_store(irec,:) = hgammar(:)
enddo
+! define and store Lagrange interpolators at all the sources
+ do i = 1,NSOURCE
+ call lagrange_any(xi_source(i),NGLLX,xigll,hxis,hpxis)
+ call lagrange_any(gamma_source(i),NGLLZ,zigll,hgammas,hpgammas)
+ hxis_store(i,:) = hxis(:)
+ hgammas_store(i,:) = hgammas(:)
+ enddo
+
! displacement, velocity, acceleration and inverse of the mass matrix for elastic elements
if(ipass == 1) then
@@ -2112,41 +2129,41 @@
allocate(b_displ_elastic(3,npoin))
allocate(b_veloc_elastic(3,npoin))
allocate(b_accel_elastic(3,npoin))
- allocate(rho_kl(npoin))
+ allocate(rho_kl(NGLLX,NGLLZ,nspec))
allocate(rho_k(npoin))
allocate(rhol_global(npoin))
- allocate(mu_kl(npoin))
+ allocate(mu_kl(NGLLX,NGLLZ,nspec))
allocate(mu_k(npoin))
allocate(mul_global(npoin))
- allocate(kappa_kl(npoin))
+ allocate(kappa_kl(NGLLX,NGLLZ,nspec))
allocate(kappa_k(npoin))
allocate(kappal_global(npoin))
- allocate(rhop_kl(npoin))
- allocate(alpha_kl(npoin))
- allocate(beta_kl(npoin))
- allocate(rhorho_el_hessian_final2(npoin))
+ allocate(rhop_kl(NGLLX,NGLLZ,nspec))
+ allocate(alpha_kl(NGLLX,NGLLZ,nspec))
+ allocate(beta_kl(NGLLX,NGLLZ,nspec))
+ allocate(rhorho_el_hessian_final2(NGLLX,NGLLZ,nspec))
allocate(rhorho_el_hessian_temp2(npoin))
- allocate(rhorho_el_hessian_final1(npoin))
+ allocate(rhorho_el_hessian_final1(NGLLX,NGLLZ,nspec))
allocate(rhorho_el_hessian_temp1(npoin))
else
allocate(b_displ_elastic(1,1))
allocate(b_veloc_elastic(1,1))
allocate(b_accel_elastic(1,1))
- allocate(rho_kl(1))
+ allocate(rho_kl(1,1,1))
allocate(rho_k(1))
allocate(rhol_global(1))
- allocate(mu_kl(1))
+ allocate(mu_kl(1,1,1))
allocate(mu_k(1))
allocate(mul_global(1))
- allocate(kappa_kl(1))
+ allocate(kappa_kl(1,1,1))
allocate(kappa_k(1))
allocate(kappal_global(1))
- allocate(rhop_kl(1))
- allocate(alpha_kl(1))
- allocate(beta_kl(1))
- allocate(rhorho_el_hessian_final2(1))
+ allocate(rhop_kl(1,1,1))
+ allocate(alpha_kl(1,1,1))
+ allocate(beta_kl(1,1,1))
+ allocate(rhorho_el_hessian_final2(1,1,1))
allocate(rhorho_el_hessian_temp2(1))
- allocate(rhorho_el_hessian_final1(1))
+ allocate(rhorho_el_hessian_final1(1,1,1))
allocate(rhorho_el_hessian_temp1(1))
endif
@@ -2178,36 +2195,36 @@
allocate(b_displw_poroelastic(NDIM,npoin))
allocate(b_velocw_poroelastic(NDIM,npoin))
allocate(b_accelw_poroelastic(NDIM,npoin))
- allocate(rhot_kl(npoin))
+ allocate(rhot_kl(NGLLX,NGLLZ,nspec))
allocate(rhot_k(npoin))
- allocate(rhof_kl(npoin))
+ allocate(rhof_kl(NGLLX,NGLLZ,nspec))
allocate(rhof_k(npoin))
- allocate(sm_kl(npoin))
+ allocate(sm_kl(NGLLX,NGLLZ,nspec))
allocate(sm_k(npoin))
- allocate(eta_kl(npoin))
+ allocate(eta_kl(NGLLX,NGLLZ,nspec))
allocate(eta_k(npoin))
- allocate(mufr_kl(npoin))
+ allocate(mufr_kl(NGLLX,NGLLZ,nspec))
allocate(mufr_k(npoin))
- allocate(B_kl(npoin))
+ allocate(B_kl(NGLLX,NGLLZ,nspec))
allocate(B_k(npoin))
- allocate(C_kl(npoin))
+ allocate(C_kl(NGLLX,NGLLZ,nspec))
allocate(C_k(npoin))
- allocate(M_kl(npoin))
+ allocate(M_kl(NGLLX,NGLLZ,nspec))
allocate(M_k(npoin))
- allocate(rhob_kl(npoin))
- allocate(rhofb_kl(npoin))
- allocate(phi_kl(npoin))
- allocate(Bb_kl(npoin))
- allocate(Cb_kl(npoin))
- allocate(Mb_kl(npoin))
- allocate(mufrb_kl(npoin))
- allocate(rhobb_kl(npoin))
- allocate(rhofbb_kl(npoin))
- allocate(phib_kl(npoin))
- allocate(cpI_kl(npoin))
- allocate(cpII_kl(npoin))
- allocate(cs_kl(npoin))
- allocate(ratio_kl(npoin))
+ allocate(rhob_kl(NGLLX,NGLLZ,nspec))
+ allocate(rhofb_kl(NGLLX,NGLLZ,nspec))
+ allocate(phi_kl(NGLLX,NGLLZ,nspec))
+ allocate(Bb_kl(NGLLX,NGLLZ,nspec))
+ allocate(Cb_kl(NGLLX,NGLLZ,nspec))
+ allocate(Mb_kl(NGLLX,NGLLZ,nspec))
+ allocate(mufrb_kl(NGLLX,NGLLZ,nspec))
+ allocate(rhobb_kl(NGLLX,NGLLZ,nspec))
+ allocate(rhofbb_kl(NGLLX,NGLLZ,nspec))
+ allocate(phib_kl(NGLLX,NGLLZ,nspec))
+ allocate(cpI_kl(NGLLX,NGLLZ,nspec))
+ allocate(cpII_kl(NGLLX,NGLLZ,nspec))
+ allocate(cs_kl(NGLLX,NGLLZ,nspec))
+ allocate(ratio_kl(NGLLX,NGLLZ,nspec))
allocate(phil_global(npoin))
allocate(mulfr_global(npoin))
allocate(etal_f_global(npoin))
@@ -2225,36 +2242,36 @@
allocate(b_displw_poroelastic(1,1))
allocate(b_velocw_poroelastic(1,1))
allocate(b_accelw_poroelastic(1,1))
- allocate(rhot_kl(1))
+ allocate(rhot_kl(1,1,1))
allocate(rhot_k(1))
- allocate(rhof_kl(1))
+ allocate(rhof_kl(1,1,1))
allocate(rhof_k(1))
- allocate(sm_kl(1))
+ allocate(sm_kl(1,1,1))
allocate(sm_k(1))
- allocate(eta_kl(1))
+ allocate(eta_kl(1,1,1))
allocate(eta_k(1))
- allocate(mufr_kl(1))
+ allocate(mufr_kl(1,1,1))
allocate(mufr_k(1))
- allocate(B_kl(1))
+ allocate(B_kl(1,1,1))
allocate(B_k(1))
- allocate(C_kl(1))
+ allocate(C_kl(1,1,1))
allocate(C_k(1))
- allocate(M_kl(1))
+ allocate(M_kl(1,1,1))
allocate(M_k(1))
- allocate(rhob_kl(1))
- allocate(rhofb_kl(1))
- allocate(phi_kl(1))
- allocate(Bb_kl(1))
- allocate(Cb_kl(1))
- allocate(Mb_kl(1))
- allocate(mufrb_kl(1))
- allocate(rhobb_kl(1))
- allocate(rhofbb_kl(1))
- allocate(phib_kl(1))
- allocate(cpI_kl(1))
- allocate(cpII_kl(1))
- allocate(cs_kl(1))
- allocate(ratio_kl(1))
+ allocate(rhob_kl(1,1,1))
+ allocate(rhofb_kl(1,1,1))
+ allocate(phi_kl(1,1,1))
+ allocate(Bb_kl(1,1,1))
+ allocate(Cb_kl(1,1,1))
+ allocate(Mb_kl(1,1,1))
+ allocate(mufrb_kl(1,1,1))
+ allocate(rhobb_kl(1,1,1))
+ allocate(rhofbb_kl(1,1,1))
+ allocate(phib_kl(1,1,1))
+ allocate(cpI_kl(1,1,1))
+ allocate(cpII_kl(1,1,1))
+ allocate(cs_kl(1,1,1))
+ allocate(ratio_kl(1,1,1))
allocate(phil_global(1))
allocate(mulfr_global(1))
allocate(etal_f_global(1))
@@ -2293,14 +2310,14 @@
allocate(b_displ_ac(2,npoin))
allocate(b_accel_ac(2,npoin))
allocate(accel_ac(2,npoin))
- allocate(rho_ac_kl(npoin))
+ allocate(rho_ac_kl(NGLLX,NGLLZ,nspec))
allocate(rhol_ac_global(npoin))
- allocate(kappa_ac_kl(npoin))
+ allocate(kappa_ac_kl(NGLLX,NGLLZ,nspec))
allocate(kappal_ac_global(npoin))
- allocate(rhop_ac_kl(npoin))
- allocate(alpha_ac_kl(npoin))
- allocate(rhorho_ac_hessian_final2(npoin))
- allocate(rhorho_ac_hessian_final1(npoin))
+ allocate(rhop_ac_kl(NGLLX,NGLLZ,nspec))
+ allocate(alpha_ac_kl(NGLLX,NGLLZ,nspec))
+ allocate(rhorho_ac_hessian_final2(NGLLX,NGLLZ,nspec))
+ allocate(rhorho_ac_hessian_final1(NGLLX,NGLLZ,nspec))
else
! allocate unused arrays with fictitious size
allocate(b_potential_acoustic(1))
@@ -2309,14 +2326,14 @@
allocate(b_displ_ac(1,1))
allocate(b_accel_ac(1,1))
allocate(accel_ac(1,1))
- allocate(rho_ac_kl(1))
+ allocate(rho_ac_kl(1,1,1))
allocate(rhol_ac_global(1))
- allocate(kappa_ac_kl(1))
+ allocate(kappa_ac_kl(1,1,1))
allocate(kappal_ac_global(1))
- allocate(rhop_ac_kl(1))
- allocate(alpha_ac_kl(1))
- allocate(rhorho_ac_hessian_final2(1))
- allocate(rhorho_ac_hessian_final1(1))
+ allocate(rhop_ac_kl(1,1,1))
+ allocate(alpha_ac_kl(1,1,1))
+ allocate(rhorho_ac_hessian_final2(1,1,1))
+ allocate(rhorho_ac_hessian_final1(1,1,1))
endif
endif
@@ -2636,12 +2653,12 @@
else
stop 'incorrect value of DISPLAY_SUBSET_OPTION'
endif
-
- call checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,ibool,kmato, &
+ call checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,permeability,ibool,kmato, &
coord,npoin,vpImin,vpImax,vpIImin,vpIImax, &
assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat,f0,t0,initialfield, &
time_function_type,coorg,xinterp,zinterp,shape2D_display,knods,simulation_title, &
- npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic,myrank,nproc,NSOURCE,poroelastic)
+ npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic,myrank,nproc,NSOURCE,poroelastic, &
+ freq0,Q0,TURN_VISCATTENUATION_ON)
! convert receiver angle to radians
anglerec = anglerec * pi / 180.d0
@@ -3291,35 +3308,6 @@
if(isolver == 2) then
if(any_elastic) then
- write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
- open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
- if(p_sv)then !P-SV waves
- do j=1,npoin
- read(55) (b_displ_elastic(i,j), i=1,NDIM), &
- (b_veloc_elastic(i,j), i=1,NDIM), &
- (b_accel_elastic(i,j), i=1,NDIM)
- enddo
- b_displ_elastic(3,:) = b_displ_elastic(2,:)
- b_displ_elastic(2,:) = ZERO
- b_veloc_elastic(3,:) = b_veloc_elastic(2,:)
- b_veloc_elastic(2,:) = ZERO
- b_accel_elastic(3,:) = b_accel_elastic(2,:)
- b_accel_elastic(2,:) = ZERO
- else !SH (membrane) waves
- do j=1,npoin
- read(55) b_displ_elastic(2,j), &
- b_veloc_elastic(2,j), &
- b_accel_elastic(2,j)
- enddo
- b_displ_elastic(1,:) = ZERO
- b_displ_elastic(3,:) = ZERO
- b_veloc_elastic(1,:) = ZERO
- b_veloc_elastic(3,:) = ZERO
- b_accel_elastic(1,:) = ZERO
- b_accel_elastic(3,:) = ZERO
- endif
- close(55)
-
write(outputname,'(a,i6.6,a)') 'snapshot_rho_kappa_mu_',myrank
open(unit = 97, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
@@ -3327,34 +3315,20 @@
open(unit = 98, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
- rho_kl(:) = ZERO
- mu_kl(:) = ZERO
- kappa_kl(:) = ZERO
+ rho_kl(:,:,:) = ZERO
+ mu_kl(:,:,:) = ZERO
+ kappa_kl(:,:,:) = ZERO
!
- rhop_kl(:) = ZERO
- beta_kl(:) = ZERO
- alpha_kl(:) = ZERO
- rhorho_el_hessian_final2(:) = ZERO
+ rhop_kl(:,:,:) = ZERO
+ beta_kl(:,:,:) = ZERO
+ alpha_kl(:,:,:) = ZERO
+ rhorho_el_hessian_final2(:,:,:) = ZERO
rhorho_el_hessian_temp2(:) = ZERO
- rhorho_el_hessian_final1(:) = ZERO
+ rhorho_el_hessian_final1(:,:,:) = ZERO
rhorho_el_hessian_temp1(:) = ZERO
endif
if(any_poroelastic) then
- write(outputname,'(a,i6.6,a)') 'lastframe_poroelastic_s',myrank,'.bin'
- open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
- write(outputname,'(a,i6.6,a)') 'lastframe_poroelastic_w',myrank,'.bin'
- open(unit=56,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
- do j=1,npoin
- read(55) (b_displs_poroelastic(i,j), i=1,NDIM), &
- (b_velocs_poroelastic(i,j), i=1,NDIM), &
- (b_accels_poroelastic(i,j), i=1,NDIM)
- read(56) (b_displw_poroelastic(i,j), i=1,NDIM), &
- (b_velocw_poroelastic(i,j), i=1,NDIM), &
- (b_accelw_poroelastic(i,j), i=1,NDIM)
- enddo
- close(55)
- close(56)
! Primary kernels
write(outputname,'(a,i6.6,a)') 'snapshot_mu_B_C_',myrank
@@ -3387,42 +3361,33 @@
open(unit = 19, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
- rhot_kl(:) = ZERO
- rhof_kl(:) = ZERO
- eta_kl(:) = ZERO
- sm_kl(:) = ZERO
- mufr_kl(:) = ZERO
- B_kl(:) = ZERO
- C_kl(:) = ZERO
- M_kl(:) = ZERO
+ rhot_kl(:,:,:) = ZERO
+ rhof_kl(:,:,:) = ZERO
+ eta_kl(:,:,:) = ZERO
+ sm_kl(:,:,:) = ZERO
+ mufr_kl(:,:,:) = ZERO
+ B_kl(:,:,:) = ZERO
+ C_kl(:,:,:) = ZERO
+ M_kl(:,:,:) = ZERO
!
- rhob_kl(:) = ZERO
- rhofb_kl(:) = ZERO
- phi_kl(:) = ZERO
- mufrb_kl(:) = ZERO
- Bb_kl(:) = ZERO
- Cb_kl(:) = ZERO
- Mb_kl(:) = ZERO
+ rhob_kl(:,:,:) = ZERO
+ rhofb_kl(:,:,:) = ZERO
+ phi_kl(:,:,:) = ZERO
+ mufrb_kl(:,:,:) = ZERO
+ Bb_kl(:,:,:) = ZERO
+ Cb_kl(:,:,:) = ZERO
+ Mb_kl(:,:,:) = ZERO
!
- rhobb_kl(:) = ZERO
- rhofbb_kl(:) = ZERO
- phib_kl(:) = ZERO
- cs_kl(:) = ZERO
- cpI_kl(:) = ZERO
- cpII_kl(:) = ZERO
- ratio_kl(:) = ZERO
+ rhobb_kl(:,:,:) = ZERO
+ rhofbb_kl(:,:,:) = ZERO
+ phib_kl(:,:,:) = ZERO
+ cs_kl(:,:,:) = ZERO
+ cpI_kl(:,:,:) = ZERO
+ cpII_kl(:,:,:) = ZERO
+ ratio_kl(:,:,:) = ZERO
endif
if(any_acoustic) then
- write(outputname,'(a,i6.6,a)') 'lastframe_acoustic',myrank,'.bin'
- open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
- do j=1,npoin
- read(55) b_potential_acoustic(j),&
- b_potential_dot_acoustic(j),&
- b_potential_dot_dot_acoustic(j)
- enddo
- close(55)
-
write(outputname,'(a,i6.6,a)') 'snapshot_rho_kappa_',myrank
open(unit = 95, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
@@ -3430,13 +3395,13 @@
open(unit = 96, file = 'OUTPUT_FILES/'//outputname,status = 'unknown',iostat=ios)
if (ios /= 0) stop 'Error writing snapshot to disk'
- rho_ac_kl(:) = ZERO
- kappa_ac_kl(:) = ZERO
+ rho_ac_kl(:,:,:) = ZERO
+ kappa_ac_kl(:,:,:) = ZERO
!
- rhop_ac_kl(:) = ZERO
- alpha_ac_kl(:) = ZERO
- rhorho_ac_hessian_final2(:) = ZERO
- rhorho_ac_hessian_final1(:) = ZERO
+ rhop_ac_kl(:,:,:) = ZERO
+ alpha_ac_kl(:,:,:) = ZERO
+ rhorho_ac_hessian_final2(:,:,:) = ZERO
+ rhorho_ac_hessian_final1(:,:,:) = ZERO
endif
endif ! if(isover == 2)
@@ -3786,10 +3751,10 @@
! Ricker (second derivative of a Gaussian) source time function
if(time_function_type(i_source) == 1) then
- source_time_function(i_source,it) = - factor(i_source) * (ONE-TWO*aval(i_source)*(time-t0(i_source))**2) * &
- exp(-aval(i_source)*(time-t0(i_source))**2)
-! source_time_function(i_source,it) = - factor(i_source) * TWO*aval(i_source)*sqrt(aval(i_source))*&
-! (time-t0(i_source))/pi * exp(-aval(i_source)*(time-t0(i_source))**2)
+! source_time_function(i_source,it) = - factor(i_source) * (ONE-TWO*aval(i_source)*(time-t0(i_source))**2) * &
+! exp(-aval(i_source)*(time-t0(i_source))**2)
+ source_time_function(i_source,it) = - factor(i_source) * TWO*aval(i_source)*sqrt(aval(i_source))*&
+ (time-t0(i_source))/pi * exp(-aval(i_source)*(time-t0(i_source))**2)
! first derivative of a Gaussian source time function
else if(time_function_type(i_source) == 2) then
@@ -4462,7 +4427,8 @@
cssquare = mul_fr/afactor
! Approximated ratio r = amplitude "w" field/amplitude "s" field (no viscous dissipation)
-! used later for kernels calculation
+! used later for wavespeed kernels calculation, which are presently implemented for inviscid case,
+! contrary to primary and density-normalized kernels, which are consistent with viscous fluid case.
gamma1 = H_biot - phil/tortl*C_biot
gamma2 = C_biot - phil/tortl*M_biot
gamma3 = phil/tortl*( M_biot*(afactor/rhol_f + phil/tortl) - C_biot)
@@ -5135,16 +5101,31 @@
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure
if(source_type(i_source) == 1) then
+
if(isolver == 1) then ! forward wavefield
- potential_dot_dot_acoustic(iglob_source(i_source)) = potential_dot_dot_acoustic(iglob_source(i_source)) &
- - source_time_function(i_source,it)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - source_time_function(i_source,it)*hlagrange
+ enddo
+ enddo
else ! backward wavefield
- b_potential_dot_dot_acoustic(iglob_source(i_source)) = b_potential_dot_dot_acoustic(iglob_source(i_source)) &
- - source_time_function(i_source,NSTEP-it+1)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ - source_time_function(i_source,NSTEP-it+1)*hlagrange
+ enddo
+ enddo
endif
+
! moment tensor
else if(source_type(i_source) == 2) then
call exit_MPI('cannot have moment tensor source in acoustic element')
+
endif
endif ! if this processor carries the source and the source element is acoustic
enddo ! do i_source=1,NSOURCE
@@ -5182,6 +5163,7 @@
b_potential_dot_acoustic = b_potential_dot_acoustic + b_deltatover2*b_potential_dot_dot_acoustic
endif
+
! free surface for an acoustic medium
if ( nelem_acoustic_surface > 0 ) then
call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
@@ -5193,9 +5175,10 @@
endif
endif
- endif
+ endif !if(any_acoustic)
+
! *********************************************************
! ************* main solver for the elastic elements
! *********************************************************
@@ -5726,25 +5709,52 @@
if(isolver == 1) then ! forward wavefield
if(p_sv) then ! P-SV calculation
- accel_elastic(1,iglob_source(i_source)) = accel_elastic(1,iglob_source(i_source)) &
- - sin(angleforce(i_source))*source_time_function(i_source,it)
- accel_elastic(3,iglob_source(i_source)) = accel_elastic(3,iglob_source(i_source)) &
- + cos(angleforce(i_source))*source_time_function(i_source,it)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) &
+ - sin(angleforce(i_source))*source_time_function(i_source,it)*hlagrange
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) &
+ + cos(angleforce(i_source))*source_time_function(i_source,it)*hlagrange
+ enddo
+ enddo
else ! SH (membrane) calculation
- accel_elastic(2,iglob_source(i_source)) = accel_elastic(2,iglob_source(i_source)) &
- + source_time_function(i_source,it)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) &
+ + source_time_function(i_source,it)*hlagrange
+ enddo
+ enddo
endif
else ! backward wavefield
if(p_sv) then ! P-SV calculation
- b_accel_elastic(1,iglob_source(i_source)) = b_accel_elastic(1,iglob_source(i_source)) &
- - sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
- b_accel_elastic(3,iglob_source(i_source)) = b_accel_elastic(3,iglob_source(i_source)) &
- + cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) &
+ - sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1) &
+ *hlagrange
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) &
+ + cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1) &
+ *hlagrange
+ enddo
+ enddo
else ! SH (membrane) calculation
- b_accel_elastic(2,iglob_source(i_source)) = b_accel_elastic(2,iglob_source(i_source)) &
- + source_time_function(i_source,NSTEP-it+1)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) &
+ + source_time_function(i_source,NSTEP-it+1)*hlagrange
+ enddo
+ enddo
+
endif
endif !endif isolver == 1
@@ -5769,7 +5779,7 @@
b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
endif
- endif
+ endif !if(any_elastic)
! ******************************************************************************************************************
@@ -5805,7 +5815,7 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
mufr_k,B_k,NSOURCE,nrec,isolver,save_forward,&
b_absorb_poro_s_left,b_absorb_poro_s_right,b_absorb_poro_s_bottom,b_absorb_poro_s_top,&
- nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax)
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,f0(1),freq0,Q0)
@@ -5828,7 +5838,7 @@
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
C_k,M_k,NSOURCE,nrec,isolver,save_forward,&
b_absorb_poro_w_left,b_absorb_poro_w_right,b_absorb_poro_w_bottom,b_absorb_poro_w_top,&
- nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax)
+ nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,f0(1),freq0,Q0)
if(save_forward .and. isolver == 1) then
@@ -6348,28 +6358,39 @@
! collocated force
if(source_type(i_source) == 1) then
if(isolver == 1) then ! forward wavefield
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
! s
- accels_poroelastic(1,iglob_source(i_source)) = accels_poroelastic(1,iglob_source(i_source)) - &
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) - hlagrange * &
(1._CUSTOM_REAL - phil/tortl)*sin(angleforce(i_source))*source_time_function(i_source,it)
- accels_poroelastic(2,iglob_source(i_source)) = accels_poroelastic(2,iglob_source(i_source)) + &
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + hlagrange * &
(1._CUSTOM_REAL - phil/tortl)*cos(angleforce(i_source))*source_time_function(i_source,it)
! w
- accelw_poroelastic(1,iglob_source(i_source)) = accelw_poroelastic(1,iglob_source(i_source)) - &
+ accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) - hlagrange * &
(1._CUSTOM_REAL - rhol_f/rhol_bar)*sin(angleforce(i_source))*source_time_function(i_source,it)
- accelw_poroelastic(2,iglob_source(i_source)) = accelw_poroelastic(2,iglob_source(i_source)) + &
+ accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + hlagrange * &
(1._CUSTOM_REAL - rhol_f/rhol_bar)*cos(angleforce(i_source))*source_time_function(i_source,it)
-
+ enddo
+ enddo
else ! backward wavefield
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ hlagrange = hxis_store(i_source,i) * hgammas_store(i_source,j)
! b_s
- b_accels_poroelastic(1,iglob_source(i_source)) = b_accels_poroelastic(1,iglob_source(i_source)) - &
+ b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) - hlagrange * &
(1._CUSTOM_REAL - phil/tortl)*sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
- b_accels_poroelastic(2,iglob_source(i_source)) = b_accels_poroelastic(2,iglob_source(i_source)) + &
+ b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + hlagrange * &
(1._CUSTOM_REAL - phil/tortl)*cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
!b_w
- b_accelw_poroelastic(1,iglob_source(i_source)) = b_accelw_poroelastic(1,iglob_source(i_source)) - &
+ b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) - hlagrange * &
(1._CUSTOM_REAL - rhol_f/rhol_bar)*sin(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
- b_accelw_poroelastic(2,iglob_source(i_source)) = b_accelw_poroelastic(2,iglob_source(i_source)) + &
+ b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) + hlagrange * &
(1._CUSTOM_REAL - rhol_f/rhol_bar)*cos(angleforce(i_source))*source_time_function(i_source,NSTEP-it+1)
+ enddo
+ enddo
endif !endif isolver == 1
endif
@@ -6396,7 +6417,7 @@
b_velocw_poroelastic = b_velocw_poroelastic + b_deltatover2*b_accelw_poroelastic
endif
- endif
+ endif !if(any_poroelastic)
!assembling the displacements on the elastic-poro boundaries
@@ -6483,7 +6504,84 @@
enddo
endif
-! kernels calculation
+! ********************************************************************************************
+! reading lastframe for adjoint/kernels calculation
+! ********************************************************************************************
+ if(it == 1 .and. isolver == 2) then
+
+! acoustic medium
+ if(any_acoustic) then
+ write(outputname,'(a,i6.6,a)') 'lastframe_acoustic',myrank,'.bin'
+ open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ do j=1,npoin
+ read(55) b_potential_acoustic(j),&
+ b_potential_dot_acoustic(j),&
+ b_potential_dot_dot_acoustic(j)
+ enddo
+ close(55)
+
+! free surface for an acoustic medium
+ if ( nelem_acoustic_surface > 0 ) then
+ call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
+ b_potential_acoustic,acoustic_surface,ibool,nelem_acoustic_surface,npoin,nspec)
+ endif
+ endif
+
+! elastic medium
+ if(any_elastic) then
+ write(outputname,'(a,i6.6,a)') 'lastframe_elastic',myrank,'.bin'
+ open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ if(p_sv)then !P-SV waves
+ do j=1,npoin
+ read(55) (b_displ_elastic(i,j), i=1,NDIM), &
+ (b_veloc_elastic(i,j), i=1,NDIM), &
+ (b_accel_elastic(i,j), i=1,NDIM)
+ enddo
+ b_displ_elastic(3,:) = b_displ_elastic(2,:)
+ b_displ_elastic(2,:) = ZERO
+ b_veloc_elastic(3,:) = b_veloc_elastic(2,:)
+ b_veloc_elastic(2,:) = ZERO
+ b_accel_elastic(3,:) = b_accel_elastic(2,:)
+ b_accel_elastic(2,:) = ZERO
+ else !SH (membrane) waves
+ do j=1,npoin
+ read(55) b_displ_elastic(2,j), &
+ b_veloc_elastic(2,j), &
+ b_accel_elastic(2,j)
+ enddo
+ b_displ_elastic(1,:) = ZERO
+ b_displ_elastic(3,:) = ZERO
+ b_veloc_elastic(1,:) = ZERO
+ b_veloc_elastic(3,:) = ZERO
+ b_accel_elastic(1,:) = ZERO
+ b_accel_elastic(3,:) = ZERO
+ endif
+ close(55)
+ endif
+
+! poroelastic medium
+ if(any_poroelastic) then
+ write(outputname,'(a,i6.6,a)') 'lastframe_poroelastic_s',myrank,'.bin'
+ open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ write(outputname,'(a,i6.6,a)') 'lastframe_poroelastic_w',myrank,'.bin'
+ open(unit=56,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ do j=1,npoin
+ read(55) (b_displs_poroelastic(i,j), i=1,NDIM), &
+ (b_velocs_poroelastic(i,j), i=1,NDIM), &
+ (b_accels_poroelastic(i,j), i=1,NDIM)
+ read(56) (b_displw_poroelastic(i,j), i=1,NDIM), &
+ (b_velocw_poroelastic(i,j), i=1,NDIM), &
+ (b_accelw_poroelastic(i,j), i=1,NDIM)
+ enddo
+ close(55)
+ close(56)
+ endif
+
+ endif ! if(it == 1 .and. isolver == 2)
+
+! ********************************************************************************************
+! kernels calculation
+! ********************************************************************************************
if(any_elastic .and. isolver == 2) then ! kernels calculation
do iglob = 1,npoin
rho_k(iglob) = accel_elastic(1,iglob)*b_displ_elastic(1,iglob) +&
@@ -6803,20 +6901,27 @@
endif
enddo
- do iglob =1,npoin
- rho_ac_kl(iglob) = rho_ac_kl(iglob) - rhol_ac_global(iglob) * &
+ do ispec = 1,nspec
+ if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ rho_ac_kl(i,j,ispec) = rho_ac_kl(i,j,ispec) - rhol_ac_global(iglob) * &
dot_product(accel_ac(:,iglob),b_displ_ac(:,iglob)) * deltat
- kappa_ac_kl(iglob) = kappa_ac_kl(iglob) - kappal_ac_global(iglob) * &
+ kappa_ac_kl(i,j,ispec) = kappa_ac_kl(i,j,ispec) - kappal_ac_global(iglob) * &
potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob) * &
b_potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob)&
* deltat
!
- rhop_ac_kl(iglob) = rho_ac_kl(iglob) + kappa_ac_kl(iglob)
- alpha_ac_kl(iglob) = TWO * kappa_ac_kl(iglob)
- rhorho_ac_hessian_final1(iglob) = rhorho_ac_hessian_final1(iglob) + &
+ rhop_ac_kl(i,j,ispec) = rho_ac_kl(i,j,ispec) + kappa_ac_kl(i,j,ispec)
+ alpha_ac_kl(i,j,ispec) = TWO * kappa_ac_kl(i,j,ispec)
+ rhorho_ac_hessian_final1(i,j,ispec) = rhorho_ac_hessian_final1(i,j,ispec) + &
dot_product(accel_ac(:,iglob),accel_ac(:,iglob)) * deltat
- rhorho_ac_hessian_final2(iglob) = rhorho_ac_hessian_final2(iglob) + &
+ rhorho_ac_hessian_final2(i,j,ispec) = rhorho_ac_hessian_final2(i,j,ispec) + &
dot_product(accel_ac(:,iglob),b_accel_ac(:,iglob)) * deltat
+ enddo
+ enddo
+ endif
enddo
endif !if(any_acoustic)
@@ -6825,40 +6930,39 @@
do ispec = 1, nspec
if(elastic(ispec)) then
- do k = 1, NGLLZ
+ do j = 1, NGLLZ
do i = 1, NGLLX
- iglob = ibool(i,k,ispec)
+ iglob = ibool(i,j,ispec)
mul_global(iglob) = poroelastcoef(2,1,kmato(ispec))
kappal_global(iglob) = poroelastcoef(3,1,kmato(ispec)) - 4._CUSTOM_REAL*mul_global(iglob)/3._CUSTOM_REAL
rhol_global(iglob) = density(1,kmato(ispec))
+
+ rho_kl(i,j,ispec) = rho_kl(i,j,ispec) - rhol_global(iglob) * rho_k(iglob) * deltat
+ mu_kl(i,j,ispec) = mu_kl(i,j,ispec) - TWO * mul_global(iglob) * mu_k(iglob) * deltat
+ kappa_kl(i,j,ispec) = kappa_kl(i,j,ispec) - kappal_global(iglob) * kappa_k(iglob) * deltat
+!
+ rhop_kl(i,j,ispec) = rho_kl(i,j,ispec) + kappa_kl(i,j,ispec) + mu_kl(i,j,ispec)
+ beta_kl(i,j,ispec) = TWO * (mu_kl(i,j,ispec) - 4._CUSTOM_REAL * mul_global(iglob) &
+ / (3._CUSTOM_REAL * kappal_global(iglob)) * kappa_kl(i,j,ispec))
+ alpha_kl(i,j,ispec) = TWO * (1._CUSTOM_REAL + 4._CUSTOM_REAL * mul_global(iglob)/&
+ (3._CUSTOM_REAL * kappal_global(iglob))) * kappa_kl(i,j,ispec)
+ rhorho_el_hessian_final1(i,j,ispec) = rhorho_el_hessian_final1(i,j,ispec) + rhorho_el_hessian_temp1(iglob) * deltat
+ rhorho_el_hessian_final2(i,j,ispec) = rhorho_el_hessian_final2(i,j,ispec) + rhorho_el_hessian_temp2(iglob) * deltat
+
enddo
enddo
endif
enddo
- do iglob =1,npoin
- rho_kl(iglob) = rho_kl(iglob) - rhol_global(iglob) * rho_k(iglob) * deltat
- mu_kl(iglob) = mu_kl(iglob) - TWO * mul_global(iglob) * mu_k(iglob) * deltat
- kappa_kl(iglob) = kappa_kl(iglob) - kappal_global(iglob) * kappa_k(iglob) * deltat
-!
- rhop_kl(iglob) = rho_kl(iglob) + kappa_kl(iglob) + mu_kl(iglob)
- beta_kl(iglob) = TWO * (mu_kl(iglob) - 4._CUSTOM_REAL * mul_global(iglob) &
- / (3._CUSTOM_REAL * kappal_global(iglob)) * kappa_kl(iglob))
- alpha_kl(iglob) = TWO * (1._CUSTOM_REAL + 4._CUSTOM_REAL * mul_global(iglob)/&
- (3._CUSTOM_REAL * kappal_global(iglob))) * kappa_kl(iglob)
- rhorho_el_hessian_final1(iglob) = rhorho_el_hessian_final1(iglob) + rhorho_el_hessian_temp1(iglob) * deltat
- rhorho_el_hessian_final2(iglob) = rhorho_el_hessian_final2(iglob) + rhorho_el_hessian_temp2(iglob) * deltat
- enddo
-
endif !if(any_elastic)
if(any_poroelastic) then
do ispec = 1, nspec
if(poroelastic(ispec)) then
- do k = 1, NGLLZ
+ do j = 1, NGLLZ
do i = 1, NGLLX
- iglob = ibool(i,k,ispec)
+ iglob = ibool(i,j,ispec)
phil_global(iglob) = porosity(kmato(ispec))
tortl_global(iglob) = tortuosity(kmato(ispec))
rhol_s_global(iglob) = density(1,kmato(ispec))
@@ -6870,62 +6974,57 @@
permlxz_global(iglob) = permeability(2,kmato(ispec))
permlzz_global(iglob) = permeability(3,kmato(ispec))
mulfr_global(iglob) = poroelastcoef(2,3,kmato(ispec))
- enddo
- enddo
- endif
- enddo
- do iglob =1,npoin
- rhot_kl(iglob) = rhot_kl(iglob) - deltat * rhol_bar_global(iglob) * rhot_k(iglob)
- rhof_kl(iglob) = rhof_kl(iglob) - deltat * rhol_f_global(iglob) * rhof_k(iglob)
- sm_kl(iglob) = sm_kl(iglob) - deltat * rhol_f_global(iglob)*tortl_global(iglob)/phil_global(iglob) * sm_k(iglob)
+ rhot_kl(i,j,ispec) = rhot_kl(i,j,ispec) - deltat * rhol_bar_global(iglob) * rhot_k(iglob)
+ rhof_kl(i,j,ispec) = rhof_kl(i,j,ispec) - deltat * rhol_f_global(iglob) * rhof_k(iglob)
+ sm_kl(i,j,ispec) = sm_kl(i,j,ispec) - deltat * rhol_f_global(iglob)*tortl_global(iglob)/phil_global(iglob) * sm_k(iglob)
!at the moment works with constant permeability
- eta_kl(iglob) = eta_kl(iglob) - deltat * etal_f_global(iglob)/permlxx_global(iglob) * eta_k(iglob)
- B_kl(iglob) = B_kl(iglob) - deltat * B_k(iglob)
- C_kl(iglob) = C_kl(iglob) - deltat * C_k(iglob)
- M_kl(iglob) = M_kl(iglob) - deltat * M_k(iglob)
- mufr_kl(iglob) = mufr_kl(iglob) - TWO * deltat * mufr_k(iglob)
+ eta_kl(i,j,ispec) = eta_kl(i,j,ispec) - deltat * etal_f_global(iglob)/permlxx_global(iglob) * eta_k(iglob)
+ B_kl(i,j,ispec) = B_kl(i,j,ispec) - deltat * B_k(iglob)
+ C_kl(i,j,ispec) = C_kl(i,j,ispec) - deltat * C_k(iglob)
+ M_kl(i,j,ispec) = M_kl(i,j,ispec) - deltat * M_k(iglob)
+ mufr_kl(i,j,ispec) = mufr_kl(i,j,ispec) - TWO * deltat * mufr_k(iglob)
! density kernels
rholb = rhol_bar_global(iglob) - phil_global(iglob)*rhol_f_global(iglob)/tortl_global(iglob)
- rhob_kl(iglob) = rhot_kl(iglob) + B_kl(iglob) + mufr_kl(iglob)
- rhofb_kl(iglob) = rhof_kl(iglob) + C_kl(iglob) + M_kl(iglob) + sm_kl(iglob)
- Bb_kl(iglob) = B_kl(iglob)
- Cb_kl(iglob) = C_kl(iglob)
- Mb_kl(iglob) = M_kl(iglob)
- mufrb_kl(iglob) = mufr_kl(iglob)
- phi_kl(iglob) = - sm_kl(iglob) - M_kl(iglob)
+ rhob_kl(i,j,ispec) = rhot_kl(i,j,ispec) + B_kl(i,j,ispec) + mufr_kl(i,j,ispec)
+ rhofb_kl(i,j,ispec) = rhof_kl(i,j,ispec) + C_kl(i,j,ispec) + M_kl(i,j,ispec) + sm_kl(i,j,ispec)
+ Bb_kl(i,j,ispec) = B_kl(i,j,ispec)
+ Cb_kl(i,j,ispec) = C_kl(i,j,ispec)
+ Mb_kl(i,j,ispec) = M_kl(i,j,ispec)
+ mufrb_kl(i,j,ispec) = mufr_kl(i,j,ispec)
+ phi_kl(i,j,ispec) = - sm_kl(i,j,ispec) - M_kl(i,j,ispec)
! wave speed kernels
dd1 = (1._CUSTOM_REAL+rholb/rhol_f_global(iglob))*ratio**2 + 2._CUSTOM_REAL*ratio +&
tortl_global(iglob)/phil_global(iglob)
- rhobb_kl(iglob) = rhob_kl(iglob) - &
+ rhobb_kl(i,j,ispec) = rhob_kl(i,j,ispec) - &
phil_global(iglob)*rhol_f_global(iglob)/(tortl_global(iglob)*B_biot) * &
(cpIIsquare + (cpIsquare - cpIIsquare)*( (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)/dd1+&
(rhol_bar_global(iglob)**2*ratio**2/rhol_f_global(iglob)**2*(phil_global(iglob)/tortl_global(iglob)*&
ratio+1)*(phil_global(iglob)/tortl_global(iglob)*ratio+phil_global(iglob)/tortl_global(iglob)*&
(1+rhol_f_global(iglob)/rhol_bar_global(iglob))-1) )/dd1**2 )- FOUR_THIRDS*cssquare )*&
- Bb_kl(iglob) - &
+ Bb_kl(i,j,ispec) - &
rhol_bar_global(iglob)*ratio**2/M_biot * (cpIsquare - cpIIsquare)* &
- (phil_global(iglob)/tortl_global(iglob)*ratio + 1._CUSTOM_REAL)**2/dd1**2*Mb_kl(iglob) + &
+ (phil_global(iglob)/tortl_global(iglob)*ratio + 1._CUSTOM_REAL)**2/dd1**2*Mb_kl(i,j,ispec) + &
rhol_bar_global(iglob)*ratio/C_biot * (cpIsquare - cpIIsquare)* (&
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)/dd1 - &
phil_global(iglob)*ratio/tortl_global(iglob)*(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*&
- (1+rhol_bar_global(iglob)*ratio/rhol_f_global(iglob))/dd1**2)*Cb_kl(iglob)+ &
- phil_global(iglob)*rhol_f_global(iglob)*cssquare/(tortl_global(iglob)*mulfr_global(iglob))*mufrb_kl(iglob)
- rhofbb_kl(iglob) = rhofb_kl(iglob) + &
+ (1+rhol_bar_global(iglob)*ratio/rhol_f_global(iglob))/dd1**2)*Cb_kl(i,j,ispec)+ &
+ phil_global(iglob)*rhol_f_global(iglob)*cssquare/(tortl_global(iglob)*mulfr_global(iglob))*mufrb_kl(i,j,ispec)
+ rhofbb_kl(i,j,ispec) = rhofb_kl(i,j,ispec) + &
phil_global(iglob)*rhol_f_global(iglob)/(tortl_global(iglob)*B_biot) * &
(cpIIsquare + (cpIsquare - cpIIsquare)*( (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)/dd1+&
(rhol_bar_global(iglob)**2*ratio**2/rhol_f_global(iglob)**2*(phil_global(iglob)/tortl_global(iglob)*&
ratio+1)*(phil_global(iglob)/tortl_global(iglob)*ratio+phil_global(iglob)/tortl_global(iglob)*&
(1+rhol_f_global(iglob)/rhol_bar_global(iglob))-1) )/dd1**2 )- FOUR_THIRDS*cssquare )*&
- Bb_kl(iglob) + &
+ Bb_kl(i,j,ispec) + &
rhol_bar_global(iglob)*ratio**2/M_biot * (cpIsquare - cpIIsquare)* &
- (phil_global(iglob)/tortl_global(iglob)*ratio + 1._CUSTOM_REAL)**2/dd1**2*Mb_kl(iglob) - &
+ (phil_global(iglob)/tortl_global(iglob)*ratio + 1._CUSTOM_REAL)**2/dd1**2*Mb_kl(i,j,ispec) - &
rhol_bar_global(iglob)*ratio/C_biot * (cpIsquare - cpIIsquare)* (&
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)/dd1 - &
phil_global(iglob)*ratio/tortl_global(iglob)*(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*&
- (1+rhol_bar_global(iglob)*ratio/rhol_f_global(iglob))/dd1**2)*Cb_kl(iglob)- &
- phil_global(iglob)*rhol_f_global(iglob)*cssquare/(tortl_global(iglob)*mulfr_global(iglob))*mufrb_kl(iglob)
- phib_kl(iglob) = phi_kl(iglob) - &
+ (1+rhol_bar_global(iglob)*ratio/rhol_f_global(iglob))/dd1**2)*Cb_kl(i,j,ispec)- &
+ phil_global(iglob)*rhol_f_global(iglob)*cssquare/(tortl_global(iglob)*mulfr_global(iglob))*mufrb_kl(i,j,ispec)
+ phib_kl(i,j,ispec) = phi_kl(i,j,ispec) - &
phil_global(iglob)*rhol_bar_global(iglob)/(tortl_global(iglob)*B_biot) * ( cpIsquare - rhol_f_global(iglob)/&
rhol_bar_global(iglob)*cpIIsquare- &
(cpIsquare-cpIIsquare)*( (TWO*ratio**2*phil_global(iglob)/tortl_global(iglob) + (1._CUSTOM_REAL+&
@@ -6934,66 +7033,70 @@
ratio/tortl_global(iglob)+phil_global(iglob)/tortl_global(iglob)*(1._CUSTOM_REAL+rhol_f_global(iglob)/&
rhol_bar_global(iglob))-1._CUSTOM_REAL)*((1._CUSTOM_REAL+rhol_bar_global(iglob)/rhol_f_global(iglob)-&
TWO*phil_global(iglob)/tortl_global(iglob))*ratio**2+TWO*ratio)/dd1**2 ) - &
- FOUR_THIRDS*rhol_f_global(iglob)*cssquare/rhol_bar_global(iglob) )*Bb_kl(iglob) + &
+ FOUR_THIRDS*rhol_f_global(iglob)*cssquare/rhol_bar_global(iglob) )*Bb_kl(i,j,ispec) + &
rhol_f_global(iglob)/M_biot * (cpIsquare-cpIIsquare)*(&
TWO*ratio*(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)/dd1 - &
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)**2*((1._CUSTOM_REAL+rhol_bar_global(iglob)/&
rhol_f_global(iglob)-TWO*phil_global(iglob)/tortl_global(iglob))*ratio**2+TWO*ratio)/dd1**2&
- )*Mb_kl(iglob) + &
+ )*Mb_kl(i,j,ispec) + &
phil_global(iglob)*rhol_f_global(iglob)/(tortl_global(iglob)*C_biot)*(cpIsquare-cpIIsquare)*ratio* (&
(1._CUSTOM_REAL+rhol_f_global(iglob)/rhol_bar_global(iglob)*ratio)/dd1 - &
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*(1._CUSTOM_REAL+rhol_bar_global(iglob)/&
rhol_f_global(iglob)*ratio)*((1._CUSTOM_REAL+rhol_bar_global(iglob)/rhol_f_global(iglob)-TWO*&
phil_global(iglob)/tortl_global(iglob))*ratio+TWO)/dd1**2&
- )*Cb_kl(iglob) -&
- phil_global(iglob)*rhol_f_global(iglob)*cssquare/(tortl_global(iglob)*mulfr_global(iglob))*mufrb_kl(iglob)
- cpI_kl(iglob) = 2._CUSTOM_REAL*cpIsquare/B_biot*rhol_bar_global(iglob)*( &
+ )*Cb_kl(i,j,ispec) -&
+ phil_global(iglob)*rhol_f_global(iglob)*cssquare/(tortl_global(iglob)*mulfr_global(iglob))*mufrb_kl(i,j,ispec)
+ cpI_kl(i,j,ispec) = 2._CUSTOM_REAL*cpIsquare/B_biot*rhol_bar_global(iglob)*( &
1._CUSTOM_REAL-phil_global(iglob)/tortl_global(iglob) + &
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*(phil_global(iglob)/tortl_global(iglob)*&
ratio+phil_global(iglob)/tortl_global(iglob)*(1._CUSTOM_REAL+rhol_f_global(iglob)/rhol_bar_global(iglob))-&
1._CUSTOM_REAL)/dd1 &
- )* Bb_kl(iglob) +&
+ )* Bb_kl(i,j,ispec) +&
2._CUSTOM_REAL*cpIsquare*rhol_f_global(iglob)*tortl_global(iglob)/(phil_global(iglob)*M_biot) *&
- (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)**2/dd1*Mb_kl(iglob)+&
+ (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)**2/dd1*Mb_kl(i,j,ispec)+&
2._CUSTOM_REAL*cpIsquare*rhol_f_global(iglob)/C_biot * &
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*(1._CUSTOM_REAL+rhol_bar_global(iglob)/&
- rhol_f_global(iglob)*ratio)/dd1*Cb_kl(iglob)
- cpII_kl(iglob) = 2._CUSTOM_REAL*cpIIsquare*rhol_bar_global(iglob)/B_biot * (&
+ rhol_f_global(iglob)*ratio)/dd1*Cb_kl(i,j,ispec)
+ cpII_kl(i,j,ispec) = 2._CUSTOM_REAL*cpIIsquare*rhol_bar_global(iglob)/B_biot * (&
phil_global(iglob)*rhol_f_global(iglob)/(tortl_global(iglob)*rhol_bar_global(iglob)) - &
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*(phil_global(iglob)/tortl_global(iglob)*&
ratio+phil_global(iglob)/tortl_global(iglob)*(1._CUSTOM_REAL+rhol_f_global(iglob)/rhol_bar_global(iglob))-&
- 1._CUSTOM_REAL)/dd1 ) * Bb_kl(iglob) +&
+ 1._CUSTOM_REAL)/dd1 ) * Bb_kl(i,j,ispec) +&
2._CUSTOM_REAL*cpIIsquare*rhol_f_global(iglob)*tortl_global(iglob)/(phil_global(iglob)*M_biot) * (&
- 1._CUSTOM_REAL - (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)**2/dd1 )*Mb_kl(iglob) + &
+ 1._CUSTOM_REAL - (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)**2/dd1 )*Mb_kl(i,j,ispec) + &
2._CUSTOM_REAL*cpIIsquare*rhol_f_global(iglob)/C_biot * (&
1._CUSTOM_REAL - (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*(1._CUSTOM_REAL+&
- rhol_bar_global(iglob)/rhol_f_global(iglob)*ratio)/dd1 )*Cb_kl(iglob)
- cs_kl(iglob) = - 8._CUSTOM_REAL/3._CUSTOM_REAL*cssquare*rhol_bar_global(iglob)/B_biot*(1._CUSTOM_REAL-&
- phil_global(iglob)*rhol_f_global(iglob)/(tortl_global(iglob)*rhol_bar_global(iglob)))*Bb_kl(iglob) + &
+ rhol_bar_global(iglob)/rhol_f_global(iglob)*ratio)/dd1 )*Cb_kl(i,j,ispec)
+ cs_kl(i,j,ispec) = - 8._CUSTOM_REAL/3._CUSTOM_REAL*cssquare*rhol_bar_global(iglob)/B_biot*(1._CUSTOM_REAL-&
+ phil_global(iglob)*rhol_f_global(iglob)/(tortl_global(iglob)*rhol_bar_global(iglob)))*Bb_kl(i,j,ispec) + &
2._CUSTOM_REAL*(rhol_bar_global(iglob)-rhol_f_global(iglob)*phil_global(iglob)/tortl_global(iglob))/&
- mulfr_global(iglob)*cssquare*mufrb_kl(iglob)
- ratio_kl(iglob) = ratio*rhol_bar_global(iglob)*phil_global(iglob)/(tortl_global(iglob)*B_biot) * &
+ mulfr_global(iglob)*cssquare*mufrb_kl(i,j,ispec)
+ ratio_kl(i,j,ispec) = ratio*rhol_bar_global(iglob)*phil_global(iglob)/(tortl_global(iglob)*B_biot) * &
(cpIsquare-cpIIsquare) * ( &
phil_global(iglob)/tortl_global(iglob)*(2._CUSTOM_REAL*ratio+1._CUSTOM_REAL+rhol_f_global(iglob)/ &
rhol_bar_global(iglob))/dd1 - (phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)*&
(phil_global(iglob)/tortl_global(iglob)*ratio+phil_global(iglob)/tortl_global(iglob)*(&
1._CUSTOM_REAL+rhol_f_global(iglob)/rhol_bar_global(iglob))-1._CUSTOM_REAL)*(2._CUSTOM_REAL*ratio*(&
1._CUSTOM_REAL+rhol_bar_global(iglob)/rhol_f_global(iglob)-phil_global(iglob)/tortl_global(iglob)) +&
- 2._CUSTOM_REAL)/dd1**2 )*Bb_kl(iglob) + &
+ 2._CUSTOM_REAL)/dd1**2 )*Bb_kl(i,j,ispec) + &
ratio*rhol_f_global(iglob)*tortl_global(iglob)/(phil_global(iglob)*M_biot)*(cpIsquare-cpIIsquare) * &
2._CUSTOM_REAL*phil_global(iglob)/tortl_global(iglob) * (&
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)/dd1 - &
(phil_global(iglob)/tortl_global(iglob)*ratio+1._CUSTOM_REAL)**2*((1._CUSTOM_REAL+rhol_bar_global(iglob)/&
- rhol_f_global(iglob)-phil_global(iglob)/tortl_global(iglob))*ratio+1._CUSTOM_REAL)/dd1**2 )*Mb_kl(iglob) +&
+ rhol_f_global(iglob)-phil_global(iglob)/tortl_global(iglob))*ratio+1._CUSTOM_REAL)/dd1**2 )*Mb_kl(i,j,ispec) +&
ratio*rhol_f_global(iglob)/C_biot*(cpIsquare-cpIIsquare) * (&
(2._CUSTOM_REAL*phil_global(iglob)*rhol_bar_global(iglob)*ratio/(tortl_global(iglob)*rhol_f_global(iglob))+&
phil_global(iglob)/tortl_global(iglob)+rhol_bar_global(iglob)/rhol_f_global(iglob))/dd1 - &
2._CUSTOM_REAL*phil_global(iglob)/tortl_global(iglob)*(phil_global(iglob)/tortl_global(iglob)*ratio+&
1._CUSTOM_REAL)*(1._CUSTOM_REAL+rhol_bar_global(iglob)/rhol_f_global(iglob)*ratio)*((1._CUSTOM_REAL+&
rhol_bar_global(iglob)/rhol_f_global(iglob)-phil_global(iglob)/tortl_global(iglob))*ratio+1._CUSTOM_REAL)/&
- dd1**2 )*Cb_kl(iglob)
- enddo
+ dd1**2 )*Cb_kl(i,j,ispec)
+ enddo
+ enddo
+ endif
+ enddo
+
endif ! if(any_poroelastic)
endif ! if(isolver == 2)
@@ -7014,43 +7117,58 @@
endif
if(any_acoustic) then
- do iglob =1,npoin
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
xx = coord(1,iglob)
zz = coord(2,iglob)
- write(95,'(5e12.3)')xx,zz,rho_ac_kl(iglob),kappa_ac_kl(iglob)
- write(96,'(5e12.3)')rhorho_ac_hessian_final1(iglob), rhorho_ac_hessian_final2(iglob),&
- rhop_ac_kl(iglob),alpha_ac_kl(iglob)
- enddo
+ write(95,'(5e11.3)')xx,zz,rho_ac_kl(i,j,ispec),kappa_ac_kl(i,j,ispec)
+ write(96,'(5e11.3)')rhorho_ac_hessian_final1(i,j,ispec), rhorho_ac_hessian_final2(i,j,ispec),&
+ rhop_ac_kl(i,j,ispec),alpha_ac_kl(i,j,ispec)
+ enddo
+ enddo
+ enddo
close(95)
close(96)
endif
if(any_elastic) then
- do iglob =1,npoin
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
xx = coord(1,iglob)
zz = coord(2,iglob)
- write(97,'(5e12.3)')xx,zz,rho_kl(iglob),kappa_kl(iglob),mu_kl(iglob)
- write(98,'(5e12.3)')rhorho_el_hessian_final1(iglob), rhorho_el_hessian_final2(iglob),&
- rhop_kl(iglob),alpha_kl(iglob),beta_kl(iglob)
- enddo
+ write(97,'(5e11.3)')xx,zz,rho_kl(i,j,ispec),kappa_kl(i,j,ispec),mu_kl(i,j,ispec)
+ write(98,'(5e11.3)')rhorho_el_hessian_final1(i,j,ispec), rhorho_el_hessian_final2(i,j,ispec),&
+ rhop_kl(i,j,ispec),alpha_kl(i,j,ispec),beta_kl(i,j,ispec)
+ enddo
+ enddo
+ enddo
close(97)
close(98)
endif
if(any_poroelastic) then
- do iglob =1,npoin
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
xx = coord(1,iglob)
zz = coord(2,iglob)
- write(144,'(5e12.3)')xx,zz,mufr_kl(iglob),B_kl(iglob),C_kl(iglob)
- write(155,'(5e12.3)')xx,zz,M_kl(iglob),rhot_kl(iglob),rhof_kl(iglob)
- write(16,'(5e12.3)')xx,zz,sm_kl(iglob),eta_kl(iglob)
- write(17,'(5e12.3)')xx,zz,mufrb_kl(iglob),Bb_kl(iglob),Cb_kl(iglob)
- write(18,'(5e12.3)')xx,zz,Mb_kl(iglob),rhob_kl(iglob),rhofb_kl(iglob)
- write(19,'(5e12.3)')xx,zz,phi_kl(iglob),eta_kl(iglob)
- write(20,'(5e12.3)')xx,zz,cpI_kl(iglob),cpII_kl(iglob),cs_kl(iglob)
- write(21,'(5e12.3)')xx,zz,rhobb_kl(iglob),rhofbb_kl(iglob),ratio_kl(iglob)
- write(22,'(5e12.3)')xx,zz,phib_kl(iglob),eta_kl(iglob)
- enddo
+ write(144,'(5e11.3)')xx,zz,mufr_kl(i,j,ispec),B_kl(i,j,ispec),C_kl(i,j,ispec)
+ write(155,'(5e11.3)')xx,zz,M_kl(i,j,ispec),rhot_kl(i,j,ispec),rhof_kl(i,j,ispec)
+ write(16,'(5e11.3)')xx,zz,sm_kl(i,j,ispec),eta_kl(i,j,ispec)
+ write(17,'(5e11.3)')xx,zz,mufrb_kl(i,j,ispec),Bb_kl(i,j,ispec),Cb_kl(i,j,ispec)
+ write(18,'(5e11.3)')xx,zz,Mb_kl(i,j,ispec),rhob_kl(i,j,ispec),rhofb_kl(i,j,ispec)
+ write(19,'(5e11.3)')xx,zz,phi_kl(i,j,ispec),eta_kl(i,j,ispec)
+ write(20,'(5e11.3)')xx,zz,cpI_kl(i,j,ispec),cpII_kl(i,j,ispec),cs_kl(i,j,ispec)
+ write(21,'(5e11.3)')xx,zz,rhobb_kl(i,j,ispec),rhofbb_kl(i,j,ispec),ratio_kl(i,j,ispec)
+ write(22,'(5e11.3)')xx,zz,phib_kl(i,j,ispec),eta_kl(i,j,ispec)
+ enddo
+ enddo
+ enddo
close(144)
close(155)
close(16)
More information about the CIG-COMMITS
mailing list