[cig-commits] r20758 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Sep 21 07:26:55 PDT 2012
Author: dkomati1
Date: 2012-09-21 07:26:55 -0700 (Fri, 21 Sep 2012)
New Revision: 20758
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
fixed a PostScript display bug
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90 2012-09-21 12:32:56 UTC (rev 20757)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90 2012-09-21 14:26:55 UTC (rev 20758)
@@ -54,7 +54,7 @@
fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges, &
fluid_poro_acoustic_ispec,fluid_poro_acoustic_iedge,num_fluid_poro_edges, &
solid_poro_poroelastic_ispec,solid_poro_poroelastic_iedge,num_solid_poro_edges, &
- myrank,nproc,ier, &
+ poroelastic,myrank,nproc,ier, &
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
@@ -96,6 +96,7 @@
integer kmato(nspec),knods(ngnod,nspec)
integer ibool(NGLLX,NGLLZ,nspec)
+ logical, dimension(nspec) :: poroelastic
double precision xinterp(pointsdisp,pointsdisp),zinterp(pointsdisp,pointsdisp)
double precision shapeint(ngnod,pointsdisp,pointsdisp)
@@ -138,7 +139,7 @@
! US letter paper or European A4
logical :: US_LETTER
- double precision convert,x1,cpIloc,xa,za,xb,zb
+ double precision convert,x1,cpIloc,xa,za,xb,zb,lambdaplus2mu,denst
double precision z1,x2,z2,d,d1,d2,dummy,theta,thetaup,thetadown
double precision :: mul_s,kappal_s,rhol_s
@@ -1662,10 +1663,19 @@
do j=1,NGLLX-subsamp_postscript,subsamp_postscript
if((vpmax-vpmin)/vpmin > 0.02d0) then
+
if(assign_external_model) then
+
x1 = (vpext(i,j,ispec)-vpmin) / (vpmax-vpmin)
+
else
+
material = kmato(ispec)
+
+ if(poroelastic(ispec)) then
+
+ ! poroelastic material
+
! get elastic parameters of current spectral element
phil = porosity(kmato(ispec))
tortl = tortuosity(kmato(ispec))
@@ -1691,8 +1701,19 @@
cfactor = phil/(tortl*rhol_f)*(H_biot*M_biot - C_biot*C_biot)
cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
cpIloc = sqrt(cpIsquare)
+
+ else
+
+ lambdaplus2mu = poroelastcoef(3,1,material)
+ denst = density(1,material)
+ cpIloc = sqrt(lambdaplus2mu/denst)
+
+ endif
+
x1 = (cpIloc-vpmin)/(vpmax-vpmin)
+
endif
+
else
x1 = 0.5d0
endif
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-21 12:32:56 UTC (rev 20757)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-09-21 14:26:55 UTC (rev 20758)
@@ -889,7 +889,9 @@
!! DK DK for add spring to stacey absorbing boundary condition
logical :: ADD_SPRING_TO_STACEY
- double precision :: x_center_spring,z_center_spring,xmin_spring,xmax_spring,zmin_spring,zmax_spring
+ double precision :: x_center_spring,z_center_spring
+ double precision :: xmin,xmax,zmin,zmax
+ double precision :: xmin_local,xmax_local,zmin_local,zmax_local
!! DK DK for horizontal periodic conditions: detect common points between left and right edges
logical :: ADD_PERIODIC_CONDITIONS
@@ -1769,11 +1771,6 @@
deallocate(perm)
-! print *
-! print *,'Xmin,Xmax of the local mesh for proc ',myrank,' = ',minval(coord(1,:)),maxval(coord(1,:))
-! print *,'Zmin,Zmax of the local mesh for proc ',myrank,' = ',minval(coord(2,:)),maxval(coord(2,:))
-! print *
-
!! DK DK for periodic conditions: detect common points between left and right edges
if(ADD_PERIODIC_CONDITIONS) then
@@ -2088,14 +2085,44 @@
endif
- xmin_spring = minval(coord(1,:))
- xmax_spring = maxval(coord(1,:))
- zmin_spring = minval(coord(2,:))
- zmax_spring = maxval(coord(2,:))
+ xmin_local = minval(coord(1,:))
+ xmax_local = maxval(coord(1,:))
+ zmin_local = minval(coord(2,:))
+ zmax_local = maxval(coord(2,:))
- x_center_spring=(xmax_spring+xmin_spring)/2.d0
- z_center_spring=(zmax_spring+zmin_spring)/2.d0
+#ifdef USE_MPI
+ call MPI_REDUCE(xmin_local, xmin, 1, MPI_DOUBLE_PRECISION, MPI_MIN, 0, MPI_COMM_WORLD, ier)
+ call MPI_REDUCE(xmax_local, xmax, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD, ier)
+ call MPI_REDUCE(zmin_local, zmin, 1, MPI_DOUBLE_PRECISION, MPI_MIN, 0, MPI_COMM_WORLD, ier)
+ call MPI_REDUCE(zmax_local, zmax, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD, ier)
+#else
+ xmin = xmin_local
+ xmax = xmax_local
+ zmin = zmin_local
+ zmax = zmax_local
+#endif
+ if(myrank == 0) then
+ write(IOUT,*)
+ write(IOUT,*) 'Xmin,Xmax of the whole mesh = ',xmin,xmax
+ write(IOUT,*) 'Zmin,Zmax of the whole mesh = ',zmin,zmax
+ write(IOUT,*)
+
+! check that no source is located outside the mesh
+ do i = 1,NSOURCES
+ if(x_source(i) < xmin) stop 'error: at least one source has x < xmin of the mesh'
+ if(x_source(i) > xmax) stop 'error: at least one source has x > xmax of the mesh'
+
+ if(z_source(i) < zmin) stop 'error: at least one source has z < zmin of the mesh'
+ if(z_source(i) > zmax) stop 'error: at least one source has z > zmax of the mesh'
+ enddo
+
+ endif
+
+! use a spring to improve the stability of the Stacey condition
+ x_center_spring = (xmax + xmin)/2.d0
+ z_center_spring = (zmax + zmin)/2.d0
+
!
!--- save the grid of points in a file
!
@@ -7996,7 +8023,7 @@
fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges,&
fluid_poro_acoustic_ispec,fluid_poro_acoustic_iedge,num_fluid_poro_edges, &
solid_poro_poroelastic_ispec,solid_poro_poroelastic_iedge,num_solid_poro_edges, &
- myrank,nproc,ier,&
+ poroelastic,myrank,nproc,ier,&
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
@@ -8040,7 +8067,7 @@
fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges,&
fluid_poro_acoustic_ispec,fluid_poro_acoustic_iedge,num_fluid_poro_edges, &
solid_poro_poroelastic_ispec,solid_poro_poroelastic_iedge,num_solid_poro_edges, &
- myrank,nproc,ier,&
+ poroelastic,myrank,nproc,ier,&
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
@@ -8084,7 +8111,7 @@
fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge,num_fluid_solid_edges, &
fluid_poro_acoustic_ispec,fluid_poro_acoustic_iedge,num_fluid_poro_edges, &
solid_poro_poroelastic_ispec,solid_poro_poroelastic_iedge,num_solid_poro_edges, &
- myrank,nproc,ier,&
+ poroelastic,myrank,nproc,ier,&
d1_coorg_send_ps_velocity_model,d2_coorg_send_ps_velocity_model, &
d1_coorg_recv_ps_velocity_model,d2_coorg_recv_ps_velocity_model, &
d1_RGB_send_ps_velocity_model,d2_RGB_send_ps_velocity_model, &
More information about the CIG-COMMITS
mailing list