[cig-commits] [commit] devel: Modif to fix problem of generate databases with DSM (79fa689)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Oct 22 17:14:33 PDT 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/0479069f470e39347345ab3d4aff453320ac453e...57f2934cebe2e3d82d758afd2da6bf77d96c901d
>---------------------------------------------------------------
commit 79fa68985ad0ae71b71df51b8f699eeded056e93
Author: Clément Durochat <c.durochat at gmail.com>
Date: Tue Oct 14 16:29:31 2014 +0200
Modif to fix problem of generate databases with DSM
>---------------------------------------------------------------
79fa68985ad0ae71b71df51b8f699eeded056e93
doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf | Bin 12601259 -> 12581411 bytes
src/generate_databases/get_model.f90 | 58 +++---
src/generate_databases/model_external_values.f90 | 234 +++++++++++++++++------
src/meshfem3D/chunk_of_earth_ReadIasp91.f90 | 22 +--
src/specfem3D/prepare_timerun.F90 | 4 +-
5 files changed, 227 insertions(+), 91 deletions(-)
diff --git a/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf b/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf
index 59de46a..c654c31 100644
Binary files a/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf and b/doc/USER_MANUAL/manual_SPECFEM3D_Cartesian.pdf differ
diff --git a/src/generate_databases/get_model.f90 b/src/generate_databases/get_model.f90
index 3b32c35..fa0483b 100644
--- a/src/generate_databases/get_model.f90
+++ b/src/generate_databases/get_model.f90
@@ -72,18 +72,9 @@
ispec_is_elastic(:) = .false.
ispec_is_poroelastic(:) = .false.
- !! WANGYI test for the benchmark of hybrid DSM-SPECFEM3D coupling
- if (COUPLE_WITH_EXTERNAL_CODE) then
- if( nundefMat_ext_mesh > 6 .or. IMODEL == IMODEL_TOMO ) then ! changed by WANGYI
- write(*,*) 'nundefMat_ext_mesh, IMODEL, IMODEL_TOMO', nundefMat_ext_mesh, IMODEL, IMODEL_TOMO ! add by WANGYI
- call model_tomography_broadcast(myrank)
- endif
-
- else
- ! prepares tomographic models if needed for elements with undefined material definitions
- if( nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO ) then
- call model_tomography_broadcast(myrank)
- endif
+ ! prepares tomographic models if needed for elements with undefined material definitions
+ if( (nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO) .and. (.not. COUPLE_WITH_EXTERNAL_CODE) ) then
+ call model_tomography_broadcast(myrank)
endif
! prepares external model values if needed
@@ -94,6 +85,23 @@
call model_salton_trough_broadcast(myrank)
end select
+!! VM VM for coupling with DSM
+!! find the # layer where the middle of the element is located
+ if (COUPLE_WITH_EXTERNAL_CODE) then
+
+ if (NGLL == 5) then
+ ! gets xyz coordinates of GLL point
+ iglob = ibool(3,3,3,1)
+ xmesh = xstore_dummy(iglob)
+ ymesh = ystore_dummy(iglob)
+ zmesh = zstore_dummy(iglob)
+ else
+ stop 'bad number of GLL points for coupling with DSM'
+ endif
+
+ call FindLayer(xmesh,ymesh,zmesh)
+ end if
+
! ! Piero, read bedrock file
! in case, see file model_interface_bedrock.f90:
! call model_bedrock_broadcast(myrank)
@@ -161,6 +169,18 @@
ymesh = ystore_dummy(iglob)
zmesh = zstore_dummy(iglob)
+ !! VM VM for coupling with DSM
+ !! find the # layer where the middle of the element is located
+ if (COUPLE_WITH_EXTERNAL_CODE) then
+
+ if (NGLL == 5) then
+ if (i==3 .and. j==3 .and. k==3) call FindLayer(xmesh,ymesh,zmesh)
+ else
+ stop 'bad number of GLL points for coupling with DSM'
+ endif
+
+ end if
+
! material index 1: associated material number
! 1 = acoustic, 2 = elastic, 3 = poroelastic, -1 = undefined tomographic
imaterial_id = mat_ext_mesh(1,ispec)
@@ -346,17 +366,9 @@
call any_all_l( ANY(ispec_is_elastic), ELASTIC_SIMULATION )
call any_all_l( ANY(ispec_is_poroelastic), POROELASTIC_SIMULATION )
- !! WANGYI test for the benchmark of hybrid DSM-SPECFEM3D coupling
- if (COUPLE_WITH_EXTERNAL_CODE) then
- if( nundefMat_ext_mesh > 6 .or. IMODEL == IMODEL_TOMO ) then ! changed by wangyi for test
- call deallocate_tomography_files()
- endif
-
- else
- ! deallocates tomographic arrays
- if( nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO ) then
- call deallocate_tomography_files()
- endif
+ ! deallocates tomographic arrays
+ if( (nundefMat_ext_mesh > 0 .or. IMODEL == IMODEL_TOMO) .and. (.not. COUPLE_WITH_EXTERNAL_CODE) ) then
+ call deallocate_tomography_files()
endif
end subroutine get_model
diff --git a/src/generate_databases/model_external_values.f90 b/src/generate_databases/model_external_values.f90
index 27b9d8d..6dac64f 100644
--- a/src/generate_databases/model_external_values.f90
+++ b/src/generate_databases/model_external_values.f90
@@ -52,6 +52,15 @@
! end type model_external_variables
! type (model_external_variables) MEXT_V
+ ! VM VM my model for DSM coupling
+ use constants
+ ! VM VM
+ double precision, dimension (:,:), allocatable :: vpv_1D,vsv_1D,density_1D
+ double precision, dimension (:), allocatable :: zlayer
+ double precision, dimension (:), allocatable :: smooth_vp,smooth_vs
+ integer :: ilayer,nlayer,ncoeff,ndeg_poly
+ double precision :: ZREF,OLON,OLAT
+
end module external_model
!
@@ -62,10 +71,12 @@
! standard routine to setup model
-! use external_model
+ use external_model
use constants
+ use generate_databases_par, only: COUPLE_WITH_EXTERNAL_CODE
+
implicit none
integer :: myrank
@@ -80,6 +91,8 @@
!
! ADD YOUR MODEL HERE
!
+ if (COUPLE_WITH_EXTERNAL_CODE) call read_external_model_coupling()
+
!---
! the variables read are declared and stored in structure MEXT_V
@@ -93,7 +106,6 @@
!
!-------------------------------------------------------------------------------------------------
!
-
!
! subroutine read_external_model()
!
@@ -110,18 +122,63 @@
!---
!
! end subroutine read_external_model
-
-
!
!-------------------------------------------------------------------------------------------------
!
+!
+ subroutine read_external_model_coupling()
+!
+ use external_model !! VM VM custom subroutine for coupling with DSM
+!
+ implicit none
+!
+!---
+!
+! ADD YOUR MODEL HERE
+!
+ character(len=256):: filename
+ integer i,cc
+ double precision aa,bb
+
+ filename = IN_DATA_FILES_PATH(1:len_trim(IN_DATA_FILES_PATH))//'/coeff_poly_deg12'
+ open(27,file=trim(filename))
+ read(27,*) ndeg_poly
+ allocate(smooth_vp(0:ndeg_poly),smooth_vs(0:ndeg_poly))
+ do i=ndeg_poly,0,-1
+ read(27,*) aa,bb,cc
+ smooth_vp(i) = aa
+ smooth_vs(i) = bb
+ !write(*,*) a,b
+ end do
+
+ filename = IN_DATA_FILES_PATH(1:len_trim(IN_DATA_FILES_PATH))//'/model_1D.in'
+ open(27,file=trim(filename))
+ read(27,*) nlayer,ncoeff
+ allocate(vpv_1D(nlayer,ncoeff))
+ allocate(vsv_1D(nlayer,ncoeff))
+ allocate(density_1D(nlayer,ncoeff))
+ allocate(zlayer(nlayer))
+ do i=1,nlayer
+ read(27,*) zlayer(i)
+ read(27,*) vpv_1D(i,:)
+ read(27,*) vsv_1D(i,:)
+ read(27,*) density_1D(i,:)
+ end do
+ read(27,*) ZREF
+ read(27,*) OLON,OLAT
+!---
+!
+ end subroutine read_external_model_coupling
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine model_external_values(xmesh,ymesh,zmesh,rho,vp,vs,qkappa_atten,qmu_atten,iflag_aniso,idomain_id )
! given a GLL point, returns super-imposed velocity model values
- use generate_databases_par,only: nspec => NSPEC_AB,ibool,HUGEVAL,TINYVAL,IDOMAIN_ELASTIC
+ use generate_databases_par,only: nspec => NSPEC_AB,ibool,HUGEVAL,TINYVAL,IDOMAIN_ELASTIC,COUPLE_WITH_EXTERNAL_CODE
use create_regions_mesh_ext_par
@@ -145,6 +202,7 @@
integer :: idomain_id
! local parameters
+ double precision :: radius
real(kind=CUSTOM_REAL) :: x,y,z
real(kind=CUSTOM_REAL) :: xmin,xmax,ymin,ymax,zmin,zmax
real(kind=CUSTOM_REAL) :: depth
@@ -161,57 +219,123 @@
y = ymesh
z = zmesh
- ! note: z coordinate will be negative below surface
- ! convention is z-axis points up
-
- ! model dimensions
- xmin = 0._CUSTOM_REAL ! minval(xstore_dummy)
- xmax = 134000._CUSTOM_REAL ! maxval(xstore_dummy)
- ymin = 0._CUSTOM_REAL !minval(ystore_dummy)
- ymax = 134000._CUSTOM_REAL ! maxval(ystore_dummy)
- zmin = 0._CUSTOM_REAL ! minval(zstore_dummy)
- zmax = 60000._CUSTOM_REAL ! maxval(zstore_dummy)
-
- ! get approximate topography elevation at target coordinates from free surface
- call get_topo_elevation_free_closest(x,y,elevation,distmin, &
- nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
- num_free_surface_faces,free_surface_ispec,free_surface_ijk)
-
-
- ! depth in Z-direction
- if( distmin < HUGEVAL ) then
- depth = elevation - z
+ if (COUPLE_WITH_EXTERNAL_CODE) then
+ call model_1D(x,y,z,rho,vp,vs,radius)
else
- depth = zmin - z
- endif
+ ! note: z coordinate will be negative below surface
+ ! convention is z-axis points up
+
+ ! model dimensions
+ xmin = 0._CUSTOM_REAL ! minval(xstore_dummy)
+ xmax = 134000._CUSTOM_REAL ! maxval(xstore_dummy)
+ ymin = 0._CUSTOM_REAL !minval(ystore_dummy)
+ ymax = 134000._CUSTOM_REAL ! maxval(ystore_dummy)
+ zmin = 0._CUSTOM_REAL ! minval(zstore_dummy)
+ zmax = 60000._CUSTOM_REAL ! maxval(zstore_dummy)
+
+ ! get approximate topography elevation at target coordinates from free surface
+ call get_topo_elevation_free_closest(x,y,elevation,distmin, &
+ nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
+ num_free_surface_faces,free_surface_ispec,free_surface_ijk)
+
+
+ ! depth in Z-direction
+ if( distmin < HUGEVAL ) then
+ depth = elevation - z
+ else
+ depth = zmin - z
+ endif
+
+ ! normalizes depth between 0 and 1
+ if( abs( zmax - zmin ) > TINYVAL ) depth = depth / (zmax - zmin)
+
+ ! initial values (in m/s and kg/m^3)
+ rho = 2691.0_CUSTOM_REAL
+ vp = 4187.5_CUSTOM_REAL
+ vs = 2151.9_CUSTOM_REAL
+
+ ! adds a velocity depth gradient
+ ! (e.g. from PREM mantle gradients:
+ ! vp : 3.9382*6371/5.5
+ ! vs : 2.3481*6371/5.5
+ ! rho : 0.6924*6371/5.5 )
+ rho = rho + 802._CUSTOM_REAL * depth
+ vp = vp + 4562._CUSTOM_REAL * depth
+ vs = vs + 2720._CUSTOM_REAL * depth
+
+ ! attenuation: PREM crust value
+ qmu_atten=600._CUSTOM_REAL
+
+ ! no Q_kappa in this model
+ qkappa_atten = 9999._CUSTOM_REAL
+
+ ! no anisotropy
+ iflag_aniso = 0
+
+ ! elastic material
+ idomain_id = IDOMAIN_ELASTIC
- ! normalizes depth between 0 and 1
- if( abs( zmax - zmin ) > TINYVAL ) depth = depth / (zmax - zmin)
-
- ! initial values (in m/s and kg/m^3)
- rho = 2691.0_CUSTOM_REAL
- vp = 4187.5_CUSTOM_REAL
- vs = 2151.9_CUSTOM_REAL
-
- ! adds a velocity depth gradient
- ! (e.g. from PREM mantle gradients:
- ! vp : 3.9382*6371/5.5
- ! vs : 2.3481*6371/5.5
- ! rho : 0.6924*6371/5.5 )
- rho = rho + 802._CUSTOM_REAL * depth
- vp = vp + 4562._CUSTOM_REAL * depth
- vs = vs + 2720._CUSTOM_REAL * depth
-
- ! attenuation: PREM crust value
- qmu_atten=600._CUSTOM_REAL
-
- ! no Q_kappa in this model
- qkappa_atten = 9999._CUSTOM_REAL
-
- ! no anisotropy
- iflag_aniso = 0
-
- ! elastic material
- idomain_id = IDOMAIN_ELASTIC
+ endif
end subroutine model_external_values
+
+!----------------------------------------------------------------
+!! !! ================= VM VM CUSTOM SUBROUTINE FOR DSM COUPLING
+!----------------------------------------------------------------
+
+ subroutine FindLayer(x,y,z)
+ use external_model
+ implicit none
+ integer il
+ double precision radius,x,y,z
+ radius = dsqrt(x**2 + y**2 + (z+zref)**2) / 1000.d0
+
+ !write(124,*) 'RADIUS ',radius,x,y,z,z+zref,zref
+ il = 1
+ do while (radius .gt. zlayer(il).and.il.lt.nlayer)
+ il = il + 1
+ end do
+ il = il - 1
+ ilayer = il
+
+ !write(124,*) 'r, i : ',z,zref,radius, ilayer
+ end subroutine FindLayer
+
+!----------------------------------------------------------------
+
+ subroutine model_1D(x_eval,y_eval,z_eval, &
+ rho_final,vp_final,vs_final,r1)
+ use external_model
+ implicit none
+ double precision r1,radius,x_eval,y_eval,z_eval
+ double precision rho,vp,vs
+ real(kind=CUSTOM_REAL) rho_final,vp_final,vs_final
+ double precision Interpol,Xtol
+
+
+ Xtol=1d-2
+
+ radius = dsqrt(x_eval**2 + y_eval**2 + (z_eval+zref)**2)
+ radius = radius / 1000.d0
+ r1=radius
+
+ ! get vp,vs and rho
+ radius = radius / zlayer(nlayer)
+ vp = Interpol(vpv_1D,ilayer,radius,nlayer)
+ vs = Interpol(vsv_1D,ilayer,radius,nlayer)
+ rho = Interpol(density_1D,ilayer,radius,nlayer)
+
+ vp_final = vp * 1000.d0
+ vs_final = vs * 1000.d0
+ rho_final = rho * 1000.d0
+
+ end subroutine model_1D
+
+!----------------------------------------------------------------
+
+ function Interpol(v,i,x,nl)
+ implicit none
+ integer i,nl
+ double precision Interpol,x,v(nl,4)
+ Interpol = v(i,1)+x*(v(i,2)+x*(v(i,3)+x*v(i,4)))
+ end function Interpol
diff --git a/src/meshfem3D/chunk_of_earth_ReadIasp91.f90 b/src/meshfem3D/chunk_of_earth_ReadIasp91.f90
index 828aa52..9827019 100644
--- a/src/meshfem3D/chunk_of_earth_ReadIasp91.f90
+++ b/src/meshfem3D/chunk_of_earth_ReadIasp91.f90
@@ -133,13 +133,13 @@ subroutine Lyfnd(r,rb,n,i)
return
end subroutine Lyfnd
-
-function Interpol(v,i,x,nl)
- implicit none
- integer i,nl
- double precision Interpol,x,v(nl,4)
- Interpol = v(i,1)+x*(v(i,2)+x*(v(i,3)+x*v(i,4)))
-end function Interpol
+!! CD CD ==> to verify if it is not necessary here
+!!! function Interpol(v,i,x,nl)
+!!! implicit none
+!!! integer i,nl
+!!! double precision Interpol,x,v(nl,4)
+!!! Interpol = v(i,1)+x*(v(i,2)+x*(v(i,3)+x*v(i,4)))
+!!! end function Interpol
function IsNewLayer(x,r,n)
implicit none
@@ -222,7 +222,7 @@ end subroutine StorePointZ
zpoint(1)=zlayer(nlayer) - Z_DEPTH_BLOCK
write(*,*) zlayer(nlayer) , Z_DEPTH_BLOCK
!! niveau de depart
- call FindLayer(ilayer,zlayer,zpoint(1),nlayer)
+ call FindLayer_in_mesh_chunk(ilayer,zlayer,zpoint(1),nlayer)
write(*,*) ' INITIALISATION calcul du niveau de depart : '
write(*,*)
write(*,*) 'zlayer : ', zlayer
@@ -326,7 +326,7 @@ end subroutine StorePointZ
ProfForGemini(ilay-1,2) = zz(ilay+1)
ProfForGemini(ilay-1,3) = 0.5d0 * (zz(ilay) + zz(ilay+1))
- call FindLayer(niveau,zlayer, ProfForGemini(ilay-1,3),nlayer)
+ call FindLayer_in_mesh_chunk(niveau,zlayer, ProfForGemini(ilay-1,3),nlayer)
Niveau_elm(ilay-1)=niveau
write(*,'(i5,2f15.3,i10)') ilay,zz(ilay),zz(ilay+1),niveau
enddo
@@ -338,7 +338,7 @@ end subroutine StorePointZ
end subroutine CalGridProf
- subroutine FindLayer(i,z,r,n)
+ subroutine FindLayer_in_mesh_chunk(i,z,r,n)
implicit none
integer i,n
double precision z(n),r
@@ -353,7 +353,7 @@ end subroutine StorePointZ
enddo
- end subroutine FindLayer
+ end subroutine FindLayer_in_mesh_chunk
diff --git a/src/specfem3D/prepare_timerun.F90 b/src/specfem3D/prepare_timerun.F90
index d3b416c..67b4882 100644
--- a/src/specfem3D/prepare_timerun.F90
+++ b/src/specfem3D/prepare_timerun.F90
@@ -276,8 +276,8 @@
if(ELASTIC_SIMULATION) then
! switches to three-component mass matrix
- !! CD CD
- if( STACEY_ABSORBING_CONDITIONS .and. .not. COUPLE_WITH_EXTERNAL_CODE ) then
+ !! CD CD !! to verify
+ if( STACEY_ABSORBING_CONDITIONS) then
! adds boundary contributions
rmassx(:) = rmass(:) + rmassx(:)
rmassy(:) = rmass(:) + rmassy(:)
More information about the CIG-COMMITS
mailing list