[cig-commits] r15559 - in seismo/2D/SPECFEM2D/trunk: . DATA
cmorency at geodynamics.org
cmorency at geodynamics.org
Wed Aug 19 10:21:52 PDT 2009
Author: cmorency
Date: 2009-08-19 10:21:51 -0700 (Wed, 19 Aug 2009)
New Revision: 15559
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
Coupled simulations acoustic/elastic/poroelastic for unstructured meshes OK
Modified: seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION 2009-08-19 00:39:42 UTC (rev 15558)
+++ seismo/2D/SPECFEM2D/trunk/DATA/CMTSOLUTION 2009-08-19 17:21:51 UTC (rev 15559)
@@ -1,13 +1,13 @@
-# source 1
-source_surf = .false. # source inside the medium or at the surface
-xs = 1600. # source location x in meters
-zs = 2900. # source location z in meters
+#source 1
+source_surf = .true. # source inside the medium or at the surface
+xs = 2000 #4000. #2000. # source location x in meters
+zs = 00 #2000. #00. # source location z in meters
source_type = 2 # elastic force or acoustic pressure = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
-f0 = 15 # dominant source frequency (Hz) if not Dirac or Heaviside
-t0 = 0.
-angleforce = 00. # angle of the source (for a force only)
+f0 = 5.0 # dominant source frequency (Hz) if not Dirac or Heaviside
+t0 = 0.0 # time shift when multi sources (if one source, must be zero)
+angleforce = 0.0 # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
-factor = 1.d10 # amplification factor
+factor = 1.d10 # amplification factor
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-08-19 00:39:42 UTC (rev 15558)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-08-19 17:21:51 UTC (rev 15559)
@@ -1,15 +1,15 @@
# title of job, and file that contains interface data
-title = Test for M2 UPPA
-interfacesfile = interfaces_reg2layers.dat
+title = 2D SEG Tests
+interfacesfile = interfaces.dat
# data concerning mesh, when generated using third-party app (more info in README)
-read_external_mesh = .false.
-mesh_file = ./DATA/Mesh_canyon/canyon_mesh_file # file containing the mesh
-nodes_coords_file = ./DATA/Mesh_canyon/canyon_nodes_coords_file # file containing the nodes coordinates
-materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
-free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
-absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
+read_external_mesh = .true.
+mesh_file = ./DATA/XSEG_small/mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/XSEG_small/nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/XSEG_small/materials_file # file containing the material number for each element
+free_surface_file = ./DATA/XSEG_small/free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/XSEG_small/absorbing_surface_file # file containing the absorbing surface
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
@@ -19,9 +19,9 @@
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
-xmax = 4800.d0 # abscissa of right side of the model
-nx = 200 # number of elements along X
-ngnod = 9 # number of control nodes per element (4 or 9)
+xmax = 16000.d0 # abscissa of right side of the model
+nx = 300 # number of elements along X
+ngnod = 4 # number of control nodes per element (4 or 9)
initialfield = .false. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
@@ -39,8 +39,8 @@
absorbleft = .true.
# time step parameters
-nt = 2600 # total number of time steps
-deltat = 0.5d-3 # duration of a time step
+nt = 30000 # total number of time steps
+deltat = 3d-4 # duration of a time step
isolver = 1 # type of simulation 1=forward 2=adjoint + kernels
# source parameters
@@ -56,32 +56,24 @@
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
# receiver line parameters for seismograms
-seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl 6=potential
-save_forward = .true. # save the last frame, needed for adjoint simulation
+seismotype = 2 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl 6=potential
+save_forward = .false. # save the last frame, needed for adjoint simulation
generate_STATIONS = .true. # creates a STATION file in ./DATA
-nreceiverlines = 2 # number of receiver lines
+nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line
-nrec = 1 # number of receivers
-xdeb = 2000. # first receiver x in meters
-zdeb = 2933.33 # first receiver z in meters
-xfin = 3777. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 1866.67 # last receiver z in meters (ignored if onlyone receiver)
+nrec = 100 # number of receivers
+xdeb = 0. # first receiver x in meters
+zdeb = 0. # first receiver z in meters
+xfin = 15000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 0. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf = .false. # receivers inside the medium or at the surface
-# second receiver line
-nrec = 1 # number of receivers
-xdeb = 2000. # first receiver x in meters
-zdeb = 1866.67 # first receiver z in meters
-xfin = 3777. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 1866.67 # last receiver z in meters (ignored if onlyone receiver)
-enreg_surf = .false. # receivers inside the medium or at the surface
-
# display parameters
-NTSTEP_BETWEEN_OUTPUT_INFO = 200 # display frequency in time steps
-output_postscript_snapshot = .true. # output Postscript snapshot of the results
+NTSTEP_BETWEEN_OUTPUT_INFO = 500 # display frequency in time steps
+output_postscript_snapshot = .false. # output Postscript snapshot of the results
output_color_image = .true. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
cutsnaps = 1. # minimum amplitude in % for snapshots
@@ -96,13 +88,33 @@
outputgrid = .false. # save the grid in a text file or not
# velocity and density models
-nbmodels = 2 # nb of different models
+nbmodels = 22 # nb of different models
# define models as I: (model_number,1,rho,Vp,Vs,0,0,Qp,Qs) or II: (model_number,2,rho,c11,c13,c33,c44,Qp,Qs) or III: (model_number,3,rhos,rhof,phi,c,kxx,kxz,kzz,Ks,Kf,Kfr,etaf,mufr,Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, & poroelastic models simultaneously
-1 3 2200.d0 1040.d0 0.4d0 2.0 1d-11 0.d0 1d-11 6.9d9 4.0d9 6.7d9 0.0d-3 3.d9 10.d0 # 2700.d0 3000.d0 1732.051d0 0 0 10.d0 10.d0 0 0 0 0 0 0
-2 1 2500.d0 3000.d0 2000.d0 0 0 10.d0 10.d0 0 0 0 0 0 0 #1558.89d0 0 0 136.d0 136.d0 0 0 0 0 0 0
+1 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+2 1 1925.8d0 2050.d0 1183.41d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+3 1 2312.d0 4500.d0 2597.58d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+4 1 1986.3d0 2200.d0 1270.05d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+5 1 1881.4d0 1950.d0 1125.67d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+6 1 1808.1d0 1800.d0 1039.03d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+7 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+8 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+9 1 1986.3d0 2200.d0 1270.06d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+10 1 1904.0d0 2000.d0 1154.55d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+11 1 1857.9d0 1900.d0 1096.80d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+12 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+13 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+14 1 1833.5d0 1850.d0 1067.91d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+15 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+16 1 2056.6d0 2400.d0 1385.52d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+17 1 2116.0d0 2600.d0 1501.10d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+18 1 2087.6d0 2500.d0 1443.35d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+19 1 1946.8d0 2100.d0 1212.27d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+20 1 1857.9d0 1900.d0 1096.80d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+21 1 2087.6d0 2500.d0 1443.35d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+22 3 2200.d0 786.3d0 0.4 2.0 1d-11 0.0 1d-11 5.341d9 2d9 3d9 0.0d-4 3.204d9 10.d0
+#22 1 1634.52d0 2272.89d0 1472.71d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 2 # nb of regions and model number for each
-1 200 1 70 2
-1 200 71 140 1
+nbregions = 1 # nb of regions and model number for each
+1 300 1 160 1
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-08-19 00:39:42 UTC (rev 15558)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-08-19 17:21:51 UTC (rev 15559)
@@ -68,6 +68,14 @@
! number=2,
! pages={368-392}}
!
+!@ARTICLE{MoTr08,
+! author={C. Morency and J. Tromp},
+! year=2008,
+! title={Spectral-element simulations of wave propagation in poroelastic media},
+! journal={Geophys. J. Int.},
+! volume=175,
+! pages={301--345}}
+!
! If you use the METIS / SCOTCH / CUBIT non-structured version, please also cite:
!
! @INPROCEEDINGS{MaKoBlLe08,
@@ -81,8 +89,28 @@
! address = {Toulouse, France},
! note = {24-27 June 2008},
! url = {http://vecpar.fe.up.pt/2008}}
-
!
+! If you use the kernel capabilities of the code, please cite
+!
+! @ARTICLE{LiTr06,
+! author={Qinya Liu and Jeroen Tromp},
+! title={Finite-frequency kernels based on adjoint methods},
+! journal={Bull. Seismol. Soc. Am.},
+! year=2006,
+! volume=96,
+! number=6,
+! pages={2383-2397},
+! doi={10.1785/0120060041}}
+!
+!@ARTICLE{MoLuTr09,
+! author={C. Morency and Y. Luo and J. Tromp},
+! year=2009,
+! title={Finite-frequency kernels for wave propagation in porous media based upon adjoint methods},
+! journal=gji,
+! doi={10.1111/j.1365-246X.2009.04332}}
+!
+!
+!
! version 5.2, Dimitri Komatitsch, Nicolas Le Goff and Roland Martin, February 2008:
! - MPI implementation of the code based on domain decomposition
! with METIS or SCOTCH
@@ -341,11 +369,14 @@
integer :: num_solid_poro_edges,num_solid_poro_edges_alloc,ispec_poroelastic,ii2,jj2
logical :: coupled_elastic_poroelastic
double precision, dimension(:,:), allocatable :: displ,veloc
- double precision :: sigma_xx,sigma_xz,sigma_zz
- double precision :: b_sigma_xx,b_sigma_xz,b_sigma_zz
+ double precision :: sigma_xx,sigma_xz,sigma_zz,sigmap
+ double precision :: b_sigma_xx,b_sigma_xz,b_sigma_zz,b_sigmap
integer, dimension(:), allocatable :: ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,&
iend_top_poro,jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro
+! for detection elastic and acoustic valences
+ integer, dimension(:), allocatable :: valence_elastic,valence_acoustic,valence_poroelastic
+
! for adjoint method
logical :: save_forward ! whether or not the last frame is saved to reconstruct the forward field
integer :: isolver ! 1 = forward wavefield, 2 = backward and adjoint wavefields and kernels
@@ -3542,10 +3573,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
@@ -4161,7 +4192,50 @@
endif
+! detecting poroelastic, elastic and acoustic global points valence
+ if(coupled_acoustic_elastic .or. coupled_acoustic_poroelastic .or. coupled_elastic_poroelastic)then
+
+ allocate(valence_elastic(npoin))
+ allocate(valence_poroelastic(npoin))
+ allocate(valence_acoustic(npoin))
+
+ valence_elastic(:) = 0
+ valence_poroelastic(:) = 0
+ valence_acoustic(:) = 0
+ do ispec = 1,nspec
+ if(elastic(ispec)) then ! the element is elastic
+ do k = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,k,ispec)
+ valence_elastic(iglob) = valence_elastic(iglob) + 1
+ enddo
+ enddo
+ elseif(poroelastic(ispec)) then ! the element is poroelastic
+ do k = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,k,ispec)
+ valence_poroelastic(iglob) = valence_poroelastic(iglob) + 1
+ enddo
+ enddo
+ else ! the element is acoustic
+ do k = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,k,ispec)
+ valence_acoustic(iglob) = valence_acoustic(iglob) + 1
+ enddo
+ enddo
+ endif
+ enddo !do ispec
+
+ else
+
+ allocate(valence_elastic(1))
+ allocate(valence_poroelastic(1))
+ allocate(valence_acoustic(1))
+
+ endif !(coupled_acoustic_elastic .or. coupled_acoustic_poroelastic .or. coupled_elastic_poroelastic)
+
#ifdef USE_MPI
if(OUTPUT_ENERGY) stop 'energy calculation only currently serial only, should add an MPI_REDUCE in parallel'
#endif
@@ -4742,14 +4816,21 @@
! compute dot product
displ_n = displ_x*nx + displ_z*nz
-! formulation with generalized potential
- weight = jacobian1D * wxgll(i)
-
+ if(valence_acoustic(iglob) == valence_elastic(iglob)) then
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + weight*displ_n
+ else
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + &
+ weight*displ_n*valence_acoustic(iglob)/2._CUSTOM_REAL
+ endif
if(isolver == 2) then
+ if(valence_acoustic(iglob) == valence_elastic(iglob)) then
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) +&
weight*(b_displ_x*nx + b_displ_z*nz)
+ else
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) + &
+ weight*(b_displ_x*nx + b_displ_z*nz)*valence_acoustic(iglob)/2._CUSTOM_REAL
+ endif
endif !if(isolver == 2) then
enddo
@@ -4842,11 +4923,22 @@
! compute dot product [u_s + w]*n
displ_n = (displ_x + displw_x)*nx + (displ_z + displw_z)*nz
+ if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + weight*displ_n
+ else
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + &
+ weight*displ_n*valence_acoustic(iglob)/2._CUSTOM_REAL
+ endif
if(isolver == 2) then
+ if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) +&
weight*((b_displ_x + b_displw_x)*nx + (b_displ_z + b_displw_z)*nz)
+ else
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) + &
+ weight*((b_displ_x + b_displw_x)*nx + (b_displ_z + b_displw_z)*nz)* &
+ valence_acoustic(iglob)/2._CUSTOM_REAL
+ endif
endif !if(isolver == 2) then
enddo
@@ -4986,15 +5078,19 @@
if(iglob /= iglob2) &
call exit_MPI( 'error in solid/porous iglob detection')
- displ(1,iglob)=(displs_poroelastic(1,iglob) &
- +displ_elastic(1,iglob2))/2.d0
- displ(2,iglob)=(displs_poroelastic(2,iglob) &
- +displ_elastic(2,iglob2))/2.d0
+ displ(1,iglob)=(valence_poroelastic(iglob)*displs_poroelastic(1,iglob)&
+ +valence_elastic(iglob)*displ_elastic(1,iglob))/ &
+ (valence_poroelastic(iglob)+valence_elastic(iglob))
+ displ(2,iglob)=(valence_poroelastic(iglob)*displs_poroelastic(2,iglob) &
+ +valence_elastic(iglob)*displ_elastic(2,iglob))/ &
+ (valence_poroelastic(iglob)+valence_elastic(iglob))
- veloc(1,iglob)=(velocs_poroelastic(1,iglob) &
- +veloc_elastic(1,iglob2))/2.d0
- veloc(2,iglob)=(velocs_poroelastic(2,iglob) &
- +veloc_elastic(2,iglob2))/2.d0
+ veloc(1,iglob)=(valence_poroelastic(iglob)*velocs_poroelastic(1,iglob) &
+ +valence_elastic(iglob)*veloc_elastic(1,iglob))/ &
+ (valence_poroelastic(iglob)+valence_elastic(iglob))
+ veloc(2,iglob)=(valence_poroelastic(iglob)*velocs_poroelastic(2,iglob) &
+ +valence_elastic(iglob)*veloc_elastic(2,iglob))/ &
+ (valence_poroelastic(iglob)+valence_elastic(iglob))
enddo
enddo
@@ -5177,12 +5273,26 @@
weight = jacobian1D * wzgll(j)
endif
+ if(valence_acoustic(iglob) == valence_elastic(iglob)) then
accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure
accel_elastic(2,iglob) = accel_elastic(2,iglob) + weight*nz*pressure
+ else
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure*&
+ valence_elastic(iglob)/2._CUSTOM_REAL
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) + weight*nz*pressure*&
+ valence_elastic(iglob)/2._CUSTOM_REAL
+ endif
if(isolver == 2) then
+ if(valence_acoustic(iglob) == valence_elastic(iglob)) then
b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure
b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) + weight*nz*b_pressure
+ else
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure*&
+ valence_elastic(iglob)/2._CUSTOM_REAL
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) + weight*nz*b_pressure*&
+ valence_elastic(iglob)/2._CUSTOM_REAL
+ endif
endif !if(isolver == 2) then
enddo
@@ -5329,16 +5439,94 @@
sigma_xz = mul_G*(duz_dxl + dux_dzl)
sigma_zz = lambdalplus2mul_G*duz_dzl + lambdal_G*dux_dxl + C_biot*(dwx_dxl + dwz_dzl)
+ sigmap = C_biot*(dux_dxl + duz_dzl) + M_biot*(dwx_dxl + dwz_dzl)
+
if(isolver == 2) then
b_sigma_xx = lambdalplus2mul_G*b_dux_dxl + lambdal_G*b_duz_dzl + C_biot*(b_dwx_dxl + b_dwz_dzl)
b_sigma_xz = mul_G*(b_duz_dxl + b_dux_dzl)
b_sigma_zz = lambdalplus2mul_G*b_duz_dzl + lambdal_G*b_dux_dxl + C_biot*(b_dwx_dxl + b_dwz_dzl)
+ b_sigmap = C_biot*(b_dux_dxl + b_duz_dzl) + M_biot*(b_dwx_dxl + b_dwz_dzl)
endif
! get point values for the elastic domain, which matches our side in the inverse direction
ii2 = ivalue(ipoin1D,iedge_elastic)
jj2 = jvalue(ipoin1D,iedge_elastic)
iglob = ibool(ii2,jj2,ispec_elastic)
+! get elastic properties
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec_elastic))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec_elastic))
+ lambdalplus2mul_relaxed = poroelastcoef(3,1,kmato(ispec_elastic))
+
+! derivative along x and along z for u_s and w
+ dux_dxi = ZERO
+ duz_dxi = ZERO
+
+ dux_dgamma = ZERO
+ duz_dgamma = ZERO
+
+ if(isolver == 2) then
+ b_dux_dxi = ZERO
+ b_duz_dxi = ZERO
+
+ b_dux_dgamma = ZERO
+ b_duz_dgamma = ZERO
+ endif
+
+! first double loop over GLL points to compute and store gradients
+! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi = dux_dxi + displ_elastic(1,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
+ duz_dxi = duz_dxi + displ_elastic(2,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
+ dux_dgamma = dux_dgamma + displ_elastic(1,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
+ duz_dgamma = duz_dgamma + displ_elastic(2,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
+
+ if(isolver == 2) then
+ b_dux_dxi = b_dux_dxi + b_displ_elastic(1,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
+ b_duz_dxi = b_duz_dxi + b_displ_elastic(2,ibool(k,jj2,ispec_elastic))*hprime_xx(ii2,k)
+ b_dux_dgamma = b_dux_dgamma + b_displ_elastic(1,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
+ b_duz_dgamma = b_duz_dgamma + b_displ_elastic(2,ibool(ii2,k,ispec_elastic))*hprime_zz(jj2,k)
+ endif
+ enddo
+
+ xixl = xix(ii2,jj2,ispec_elastic)
+ xizl = xiz(ii2,jj2,ispec_elastic)
+ gammaxl = gammax(ii2,jj2,ispec_elastic)
+ gammazl = gammaz(ii2,jj2,ispec_elastic)
+
+! derivatives of displacement
+ dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+
+ duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
+ duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
+
+ if(isolver == 2) then
+ b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
+ b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
+
+ b_duz_dxl = b_duz_dxi*xixl + b_duz_dgamma*gammaxl
+ b_duz_dzl = b_duz_dxi*xizl + b_duz_dgamma*gammazl
+ endif
+! compute stress tensor
+! full anisotropy
+ if(TURN_ANISOTROPY_ON) then
+! implement anisotropy in 2D
+ sigma_xx = sigma_xx + c11val*dux_dxl + c15val*(duz_dxl + dux_dzl) + c13val*duz_dzl
+ sigma_zz = sigma_zz + c13val*dux_dxl + c35val*(duz_dxl + dux_dzl) + c33val*duz_dzl
+ sigma_xz = sigma_xz + c15val*dux_dxl + c55val*(duz_dxl + dux_dzl) + c35val*duz_dzl
+ else
+! no attenuation
+ sigma_xx = sigma_xx + lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
+ sigma_xz = sigma_xz + mul_relaxed*(duz_dxl + dux_dzl)
+ sigma_zz = sigma_zz + lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
+ endif
+
+ if(isolver == 2) then
+ b_sigma_xx = b_sigma_xx + lambdalplus2mul_relaxed*b_dux_dxl + lambdal_relaxed*b_duz_dzl
+ b_sigma_xz = b_sigma_xz + mul_relaxed*(b_duz_dxl + b_dux_dzl)
+ b_sigma_zz = b_sigma_zz + lambdalplus2mul_relaxed*b_duz_dzl + lambdal_relaxed*b_dux_dxl
+ endif
+
! compute the 1D Jacobian and the normal to the edge: for their expression see for instance
! O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics,
! Sixth Edition, electronic version, www.amazon.com, p. 204 and Figure 7.7(a),
@@ -5374,18 +5562,34 @@
weight = jacobian1D * wzgll(j)
endif
+ if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
accel_elastic(1,iglob) = accel_elastic(1,iglob) - weight* &
- (sigma_xx*nx + sigma_xz*nz)
+ (sigma_xx*nx + sigma_xz*nz +phil*sigmap*nx)/3.d0
accel_elastic(2,iglob) = accel_elastic(2,iglob) - weight* &
- (sigma_xz*nx + sigma_zz*nz)
+ (sigma_xz*nx + sigma_zz*nz +phil*sigmap*nz)/3.d0
+ else
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - weight* &
+ (sigma_xx*nx + sigma_xz*nz +phil*sigmap*nx)/3.d0/2.d0*valence_elastic(iglob)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - weight* &
+ (sigma_xz*nx + sigma_zz*nz +phil*sigmap*nz)/3.d0/2.d0*valence_elastic(iglob)
+ endif
+
if(isolver == 2) then
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - weight*( &
- b_sigma_xx*nx + b_sigma_xz*nz)
+ if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - weight* &
+ (b_sigma_xx*nx + b_sigma_xz*nz +phil*b_sigmap*nx)/3.d0
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - weight*( &
- b_sigma_xz*nx + b_sigma_zz*nz)
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - weight* &
+ (b_sigma_xz*nx + b_sigma_zz*nz +phil*b_sigmap*nz)/3.d0
+ else
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - weight* &
+ (b_sigma_xx*nx + b_sigma_xz*nz +phil*b_sigmap*nx)/3.d0/2.d0*valence_elastic(iglob)
+
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - weight* &
+ (b_sigma_xz*nx + b_sigma_zz*nz +phil*b_sigmap*nz)/3.d0/2.d0*valence_elastic(iglob)
+ endif
endif !if(isolver == 2) then
enddo
@@ -5668,6 +5872,7 @@
weight = jacobian1D * wzgll(j)
endif
+ if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
! contribution to the solid phase
accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-phil/tortl)
accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-phil/tortl)
@@ -5675,8 +5880,22 @@
! contribution to the fluid phase
accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
+ else
+! contribution to the solid phase
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-phil/tortl)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-phil/tortl)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+! contribution to the fluid phase
+ accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) + weight*nx*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+ accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + weight*nz*pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+ endif
+
if(isolver == 2) then
+ if(valence_acoustic(iglob) == valence_poroelastic(iglob)) then
! contribution to the solid phase
b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-phil/tortl)
b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-phil/tortl)
@@ -5684,6 +5903,19 @@
! contribution to the fluid phase
b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)
+ else
+! contribution to the solid phase
+ b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-phil/tortl)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+ b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-phil/tortl)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+
+! contribution to the fluid phase
+ b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) + weight*nx*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+ b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) + weight*nz*b_pressure*(1._CUSTOM_REAL-rhol_f/rhol_bar)*&
+ valence_poroelastic(iglob)/2._CUSTOM_REAL
+ endif
endif !if(isolver == 2) then
enddo ! do ipoin1D = 1,NGLLX
@@ -5717,14 +5949,6 @@
j = jvalue_inverse(ipoin1D,iedge_elastic)
iglob = ibool(i,j,ispec_elastic)
-! get poroelastic medium properties
- phil = porosity(kmato(ispec_poroelastic))
- tortl = tortuosity(kmato(ispec_poroelastic))
-!
- rhol_s = density(1,kmato(ispec_poroelastic))
- rhol_f = density(2,kmato(ispec_poroelastic))
- rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
-
! get elastic properties
lambdal_relaxed = poroelastcoef(1,1,kmato(ispec_elastic))
mul_relaxed = poroelastcoef(2,1,kmato(ispec_elastic))
@@ -5806,6 +6030,129 @@
j = jvalue(ipoin1D,iedge_poroelastic)
iglob = ibool(i,j,ispec_poroelastic)
+! get poroelastic domain paramters
+ phil = porosity(kmato(ispec_poroelastic))
+ tortl = tortuosity(kmato(ispec_poroelastic))
+!solid properties
+ mul_s = poroelastcoef(2,1,kmato(ispec_poroelastic))
+ kappal_s = poroelastcoef(3,1,kmato(ispec_poroelastic)) - 4._CUSTOM_REAL*mul_s/3._CUSTOM_REAL
+ rhol_s = density(1,kmato(ispec_poroelastic))
+!fluid properties
+ kappal_f = poroelastcoef(1,2,kmato(ispec_poroelastic))
+ rhol_f = density(2,kmato(ispec_poroelastic))
+!frame properties
+ mul_fr = poroelastcoef(2,3,kmato(ispec_poroelastic))
+ kappal_fr = poroelastcoef(3,3,kmato(ispec_poroelastic)) - 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+!Biot coefficients for the input phi
+ D_biot = kappal_s*(1._CUSTOM_REAL + phil*(kappal_s/kappal_f - 1._CUSTOM_REAL))
+ 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)
+ mul_G = mul_fr
+ lambdal_G = H_biot - 2._CUSTOM_REAL*mul_fr
+ lambdalplus2mul_G = lambdal_G + TWO*mul_G
+
+! derivative along x and along z for u_s and w
+ dux_dxi = ZERO
+ duz_dxi = ZERO
+
+ dux_dgamma = ZERO
+ duz_dgamma = ZERO
+
+ dwx_dxi = ZERO
+ dwz_dxi = ZERO
+
+ dwx_dgamma = ZERO
+ dwz_dgamma = ZERO
+
+ if(isolver == 2) then
+ b_dux_dxi = ZERO
+ b_duz_dxi = ZERO
+
+ b_dux_dgamma = ZERO
+ b_duz_dgamma = ZERO
+
+ b_dwx_dxi = ZERO
+ b_dwz_dxi = ZERO
+
+ b_dwx_dgamma = ZERO
+ b_dwz_dgamma = ZERO
+ endif
+
+! first double loop over GLL points to compute and store gradients
+! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi = dux_dxi + displs_poroelastic(1,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ duz_dxi = duz_dxi + displs_poroelastic(2,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ dux_dgamma = dux_dgamma + displs_poroelastic(1,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displs_poroelastic(2,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+
+ dwx_dxi = dwx_dxi + displw_poroelastic(1,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ dwz_dxi = dwz_dxi + displw_poroelastic(2,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ dwx_dgamma = dwx_dgamma + displw_poroelastic(1,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+ dwz_dgamma = dwz_dgamma + displw_poroelastic(2,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+ if(isolver == 2) then
+ b_dux_dxi = b_dux_dxi + b_displs_poroelastic(1,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ b_duz_dxi = b_duz_dxi + b_displs_poroelastic(2,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ b_dux_dgamma = b_dux_dgamma + b_displs_poroelastic(1,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+ b_duz_dgamma = b_duz_dgamma + b_displs_poroelastic(2,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+
+ b_dwx_dxi = b_dwx_dxi + b_displw_poroelastic(1,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ b_dwz_dxi = b_dwz_dxi + b_displw_poroelastic(2,ibool(k,j,ispec_poroelastic))*hprime_xx(i,k)
+ b_dwx_dgamma = b_dwx_dgamma + b_displw_poroelastic(1,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+ b_dwz_dgamma = b_dwz_dgamma + b_displw_poroelastic(2,ibool(i,k,ispec_poroelastic))*hprime_zz(j,k)
+ endif
+ enddo
+
+ xixl = xix(i,j,ispec_poroelastic)
+ xizl = xiz(i,j,ispec_poroelastic)
+ gammaxl = gammax(i,j,ispec_poroelastic)
+ gammazl = gammaz(i,j,ispec_poroelastic)
+
+! derivatives of displacement
+ dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+
+ duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
+ duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
+
+ dwx_dxl = dwx_dxi*xixl + dwx_dgamma*gammaxl
+ dwx_dzl = dwx_dxi*xizl + dwx_dgamma*gammazl
+
+ dwz_dxl = dwz_dxi*xixl + dwz_dgamma*gammaxl
+ dwz_dzl = dwz_dxi*xizl + dwz_dgamma*gammazl
+
+ if(isolver == 2) then
+ b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
+ b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
+
+ b_duz_dxl = b_duz_dxi*xixl + b_duz_dgamma*gammaxl
+ b_duz_dzl = b_duz_dxi*xizl + b_duz_dgamma*gammazl
+
+ b_dwx_dxl = b_dwx_dxi*xixl + b_dwx_dgamma*gammaxl
+ b_dwx_dzl = b_dwx_dxi*xizl + b_dwx_dgamma*gammazl
+
+ b_dwz_dxl = b_dwz_dxi*xixl + b_dwz_dgamma*gammaxl
+ b_dwz_dzl = b_dwz_dxi*xizl + b_dwz_dgamma*gammazl
+ endif
+! compute stress tensor
+
+! no attenuation
+ sigma_xx = sigma_xx + lambdalplus2mul_G*dux_dxl + lambdal_G*duz_dzl + C_biot*(dwx_dxl + dwz_dzl)
+ sigma_xz = sigma_xz + mul_G*(duz_dxl + dux_dzl)
+ sigma_zz = sigma_zz + lambdalplus2mul_G*duz_dzl + lambdal_G*dux_dxl + C_biot*(dwx_dxl + dwz_dzl)
+
+ sigmap = C_biot*(dux_dxl + duz_dzl) + M_biot*(dwx_dxl + dwz_dzl)
+
+ if(isolver == 2) then
+ b_sigma_xx = b_sigma_xx + lambdalplus2mul_G*b_dux_dxl + lambdal_G*b_duz_dzl + C_biot*(b_dwx_dxl + b_dwz_dzl)
+ b_sigma_xz = b_sigma_xz + mul_G*(b_duz_dxl + b_dux_dzl)
+ b_sigma_zz = b_sigma_zz + lambdalplus2mul_G*b_duz_dzl + lambdal_G*b_dux_dxl + C_biot*(b_dwx_dxl + b_dwz_dzl)
+ b_sigmap = C_biot*(b_dux_dxl + b_duz_dzl) + M_biot*(b_dwx_dxl + b_dwz_dzl)
+ endif
+
! compute the 1D Jacobian and the normal to the edge: for their expression see for instance
! O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics,
! Sixth Edition, electronic version, www.amazon.com, p. 204 and Figure 7.7(a),
@@ -5841,34 +6188,54 @@
weight = jacobian1D * wzgll(j)
endif
+ if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
! contribution to the solid phase
accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + &
- weight*(sigma_xx*nx + sigma_xz*nz)*(1._CUSTOM_REAL - phil/tortl )
+ weight*((sigma_xx+phil*sigmap)*nx + sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)
accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + &
- weight*(sigma_xz*nx + sigma_zz*nz)*(1._CUSTOM_REAL - phil/tortl )
+ weight*(sigma_xz*nx + (sigma_zz+phil*sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)
! contribution to the fluid phase
- accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) - &
- weight*(rhol_f/rhol_bar - 1._CUSTOM_REAL)*(sigma_xx*nx+sigma_xz*nz)
+! w = 0
+ else
+! contribution to the solid phase
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + &
+ weight*((sigma_xx+phil*sigmap)*nx + sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)&
+ /2.d0*valence_poroelastic(iglob)
- accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) - &
- weight*(rhol_f/rhol_bar - 1._CUSTOM_REAL)*(sigma_xz*nx+sigma_zz*nz)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + &
+ weight*(sigma_xz*nx + (sigma_zz+phil*sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)&
+ /2.d0*valence_poroelastic(iglob)
+! contribution to the fluid phase
+! w = 0
+ endif
+
if(isolver == 2) then
+ if(valence_poroelastic(iglob) == valence_elastic(iglob)) then
! contribution to the solid phase
b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + &
- weight*(b_sigma_xx*nx + b_sigma_xz*nz)*(1._CUSTOM_REAL - phil/tortl)
+ weight*((b_sigma_xx+phil*b_sigmap)*nx + b_sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)
b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + &
- weight*(b_sigma_xz*nx + b_sigma_zz*nz)*(1._CUSTOM_REAL - phil/tortl)
+ weight*(b_sigma_xz*nx + (b_sigma_zz+phil*b_sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)
! contribution to the fluid phase
- b_accelw_poroelastic(1,iglob) = b_accelw_poroelastic(1,iglob) - &
- weight*(rhol_f/rhol_bar - 1._CUSTOM_REAL)*(b_sigma_xx*nx + b_sigma_xz*nz)
+! w = 0
+ else
+! contribution to the solid phase
+ b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + &
+ weight*((b_sigma_xx+phil*b_sigmap)*nx + b_sigma_xz*nz)/3.d0*(1.d0 -1.d0/tortl)&
+ /2.d0*valence_poroelastic(iglob)
- b_accelw_poroelastic(2,iglob) = b_accelw_poroelastic(2,iglob) - &
- weight*(rhol_f/rhol_bar - 1._CUSTOM_REAL)*(b_sigma_xz*nx + b_sigma_zz*nz)
+ b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + &
+ weight*(b_sigma_xz*nx + (b_sigma_zz+phil*b_sigmap)*nz)/3.d0*(1.d0 -1.d0/tortl)&
+ /2.d0*valence_poroelastic(iglob)
+
+! contribution to the fluid phase
+! w = 0
+ endif
endif !if(isolver == 2) then
enddo
More information about the CIG-COMMITS
mailing list