[cig-commits] [commit] devel: Fixes sep model after testing it. (db88fd8)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Oct 2 14:07:09 PDT 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/858a86e5e8e4c076696d24c8c9942e93d3d663cb...7ef78c5b45b052afca446255945d0c62788f8e93
>---------------------------------------------------------------
commit db88fd82b4dfdc498a035061bf741a1d790e932d
Author: Matthieu Lefebvre <ml15 at princeton.edu>
Date: Wed Oct 1 12:02:09 2014 -0400
Fixes sep model after testing it.
>---------------------------------------------------------------
db88fd82b4dfdc498a035061bf741a1d790e932d
src/auxiliaries/combine_vol_data.f90 | 5 ++-
src/generate_databases/get_model.f90 | 4 +-
src/generate_databases/model_sep.f90 | 74 +++++++++++++++++++-----------------
src/generate_databases/parse_sep.c | 10 ++---
4 files changed, 50 insertions(+), 43 deletions(-)
diff --git a/src/auxiliaries/combine_vol_data.f90 b/src/auxiliaries/combine_vol_data.f90
index 011761d..31d3d70 100644
--- a/src/auxiliaries/combine_vol_data.f90
+++ b/src/auxiliaries/combine_vol_data.f90
@@ -104,7 +104,7 @@
logical :: STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE
logical :: ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical :: PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID,COUPLE_WITH_DSM
- character(len=MAX_STRING_LEN) :: LOCAL_PATH,TOMOGRAPHY_PATH,TRACTION_PATH
+ character(len=MAX_STRING_LEN) :: LOCAL_PATH,TOMOGRAPHY_PATH,TRACTION_PATH, SEP_MODEL_DIRECTORY
integer :: IMODEL
! ADIOS parameters
@@ -166,7 +166,8 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRACTION_PATH,COUPLE_WITH_DSM)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL, SEP_MODEL_DIRECTORY, &
+ FULL_ATTENUATION_SOLID,TRACTION_PATH,COUPLE_WITH_DSM)
if (ADIOS_FOR_MESH) then
call init_adios(value_file_name, mesh_file_name, value_handle, mesh_handle)
diff --git a/src/generate_databases/get_model.f90 b/src/generate_databases/get_model.f90
index 17d8875..ef61b23 100644
--- a/src/generate_databases/get_model.f90
+++ b/src/generate_databases/get_model.f90
@@ -383,7 +383,7 @@
use generate_databases_par,only: IMODEL, &
IMODEL_DEFAULT,IMODEL_GLL,IMODEL_1D_PREM,IMODEL_1D_CASCADIA,IMODEL_1D_SOCAL, &
IMODEL_SALTON_TROUGH,IMODEL_TOMO,IMODEL_USER_EXTERNAL,IMODEL_IPATI,IMODEL_IPATI_WATER, &
- IMODEL_1D_PREM_PB,IMODEL_GLL, &
+ IMODEL_1D_PREM_PB,IMODEL_GLL, IMODEL_SEP, &
IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,ATTENUATION_COMP_MAXIMUM
use create_regions_mesh_ext_par
@@ -426,7 +426,7 @@
! selects chosen velocity model
select case( IMODEL )
- case( IMODEL_DEFAULT,IMODEL_GLL,IMODEL_IPATI,IMODEL_IPATI_WATER )
+ case( IMODEL_DEFAULT,IMODEL_GLL,IMODEL_IPATI,IMODEL_IPATI_WATER, IMODEL_SEP )
! material values determined by mesh properties
call model_default(materials_ext_mesh,nmat_ext_mesh, &
undef_mat_prop,nundefMat_ext_mesh, &
diff --git a/src/generate_databases/model_sep.f90 b/src/generate_databases/model_sep.f90
index 8df9405..3822aa1 100644
--- a/src/generate_databases/model_sep.f90
+++ b/src/generate_databases/model_sep.f90
@@ -1,25 +1,16 @@
module model_sep_mod
implicit none
- !interface
- !subroutine parse_sep_header(header_name, &
- !nx, ny, nz, &
- !ox, oy, oz, &
- !dx, dy, dz, &
- !bin_name)
- !character(len=:) :: header_name, bin_name
- !integer :: nx, ny, nz
- !real :: ox, oy, ox, dx, dy, dz
- !end subroutine parse_sep_header
- !end interface
contains
!==============================================================================
subroutine model_sep()
use generate_databases_par, only: NGLLX, NGLLY, NGLLZ, NSPEC=>NSPEC_AB, &
- SEP_MODEL_DIRECTORY
+ SEP_MODEL_DIRECTORY, FOUR_THIRDS
use create_regions_mesh_ext_par, only: rhostore, rho_vp, rho_vs, &
- ispec_is_acoustic, ispec_is_elastic
+ ispec_is_acoustic, ispec_is_elastic, &
+ kappastore, mustore, &
+ rho_vpI, rho_vsI, rhoarraystore
real(kind=4), allocatable, dimension(:,:,:) :: vp_sep, vs_sep, rho_sep
integer :: NX, NY, NZ
@@ -34,9 +25,12 @@ subroutine model_sep()
logical :: vp_exists, vs_exists, rho_exists
! Get files from SEP_MODEL_DIRECTORY
- write(sep_header_name_vp, *) trim(SEP_MODEL_DIRECTORY) // "/vp.H"
- write(sep_header_name_vs, *) trim(SEP_MODEL_DIRECTORY) // "/vs.H"
- write(sep_header_name_rho, *) trim(SEP_MODEL_DIRECTORY) // "/rho.H"
+ sep_header_name_vp = trim(SEP_MODEL_DIRECTORY) // "/vp.H"
+ sep_header_name_vs = trim(SEP_MODEL_DIRECTORY) // "/vs.H"
+ sep_header_name_rho = trim(SEP_MODEL_DIRECTORY) // "/rho.H"
+ !write(sep_header_name_vp, '(a)') trim(SEP_MODEL_DIRECTORY) // "/vp.H"
+ !write(sep_header_name_vs, *) trim(SEP_MODEL_DIRECTORY) // "/vs.H"
+ !write(sep_header_name_rho, *) trim(SEP_MODEL_DIRECTORY) // "/rho.H"
inquire(file=trim(sep_header_name_vp), exist=vp_exists)
if (.not. vp_exists) stop "SEP vp model should exist"
@@ -53,14 +47,14 @@ subroutine model_sep()
! -> Check for unit* to be meter
! Parse only one of the header, vp is the most likely to be present
! It might be useful to make sure that all SEP files have coherent dimensions.
- call parse_sep_header(trim(sep_header_name_vp), &
- NX, NY, NZ, OX, OY, OZ, DX, DY, DZ, &
+ call parse_sep_header(trim(sep_header_name_vp) // char(0), &
+ NX, NY, NZ, OX, OY, OZ, DX, DY, DZ, &
sep_bin_vp)
if (vs_exists) then
- call parse_sep_header(trim(sep_header_name_vs), &
- NX_alt, NY_alt, NZ_alt, &
- OX_alt, OY_alt, OZ_alt, &
- DX_alt, DY_alt, DZ_alt, &
+ call parse_sep_header(trim(sep_header_name_vs) // char(0), &
+ NX_alt, NY_alt, NZ_alt, &
+ OX_alt, OY_alt, OZ_alt, &
+ DX_alt, DY_alt, DZ_alt, &
sep_bin_vs)
if ((NX /= NX_alt) .and. (NY /= NX_alt) .and. (NZ /= NZ_alt) .and. &
(OX /= OX_alt) .and. (OY /= OX_alt) .and. (OZ /= OZ_alt) .and. &
@@ -69,10 +63,10 @@ subroutine model_sep()
endif
endif
if (rho_exists) then
- call parse_sep_header(trim(sep_header_name_rho), &
- NX_alt, NY_alt, NZ_alt, &
- OX_alt, OY_alt, OZ_alt, &
- DX_alt, DY_alt, DZ_alt, &
+ call parse_sep_header(trim(sep_header_name_rho) // char(0), &
+ NX_alt, NY_alt, NZ_alt, &
+ OX_alt, OY_alt, OZ_alt, &
+ DX_alt, DY_alt, DZ_alt, &
sep_bin_rho)
if ((NX /= NX_alt) .and. (NY/= NX_alt) .and. (NZ/= NZ_alt) .and. &
(OX /= OX_alt) .and. (OY/= OX_alt) .and. (OZ/= OZ_alt) .and. &
@@ -93,6 +87,7 @@ subroutine model_sep()
call read_sep_binary_mpiio(sep_bin_vp, NX, NY, NZ, ni, nj, NZ, &
imin, jmin, kmin, vp_sep)
! Interpolate SEP values on meshfem mesh.
+ rho_vp = 0.0
call interpolate_sep_on_mesh(vp_sep, xmin, ymin, ni, nj, NZ, &
DX, DY, DZ, rho_vp)
@@ -101,10 +96,10 @@ subroutine model_sep()
!---------'
if (vs_exists) then
allocate(vs_sep(ni, nj, NZ))
- call read_sep_binary_mpiio(sep_bin_vs, NX, NY, NZ, ni, nj, &
- imin, jmin, kmin, NZ, vs_sep)
- call interpolate_sep_on_mesh(rho_sep, xmin, ymin, ni, nj, NZ, &
- DX, DY, DZ, rhostore)
+ call read_sep_binary_mpiio(sep_bin_vs, NX, NY, NZ, ni, nj, NZ, &
+ imin, jmin, kmin, vs_sep)
+ call interpolate_sep_on_mesh(vs_sep, xmin, ymin, ni, nj, NZ, &
+ DX, DY, DZ, rho_vs)
endif
!----------.
@@ -112,10 +107,10 @@ subroutine model_sep()
!----------'
if (rho_exists) then
allocate(rho_sep(ni, nj, NZ))
- call read_sep_binary_mpiio(sep_bin_rho, NX, NY, NZ, ni, nj, &
- imin, jmin, kmin, NZ, rho_sep)
- call interpolate_sep_on_mesh(vs_sep, xmin, ymin, ni, nj, NZ, &
- DX, DY, DZ, rho_vs)
+ call read_sep_binary_mpiio(sep_bin_rho, NX, NY, NZ, ni, nj, NZ, &
+ imin, jmin, kmin, rho_sep)
+ call interpolate_sep_on_mesh(rho_sep, xmin, ymin, ni, nj, NZ, &
+ DX, DY, DZ, rhostore)
endif
! Check that values around the acoustic / elastic interface is correct.
@@ -125,11 +120,20 @@ subroutine model_sep()
call correct_sep_interface()
endif ! vs_exists
+ kappastore = rhostore *( rho_vp * rho_vp - FOUR_THIRDS*rho_vs*rho_vs )
+ mustore = rhostore*rho_vs*rho_vs
+
! SPECFEM expects rho*velocity
rho_vp = rho_vp * rhostore
if (vs_exists) then
rho_vs = rho_vs * rhostore
endif
+ ! Stacey, a completer par la suite
+ !
+ rho_vpI = rho_vp
+ rho_vsI = rho_vs
+ rhoarraystore(1, :, :, :, :) = rhostore
+ rhoarraystore(2, :, :, :, :) = rhostore
if (allocated(vp_sep)) deallocate(vp_sep)
if (allocated(vs_sep)) deallocate(vs_sep)
@@ -169,6 +173,8 @@ subroutine read_sep_binary_mpiio(filename, NX, NY, NZ, ni, nj, nk, &
local_sizes = (/ ni, nj, nk /)
starting_idx = (/ imin-1, jmin-1, kmin-1 /)
+ call MPI_BARRIER(MPI_COMM_WORLD, ier)
+
! Create the 3D stencil for values consecutive in X but disjoint in Y and Z.
call MPI_Type_create_subarray(3, global_sizes, local_sizes, starting_idx, &
MPI_ORDER_FORTRAN, MPI_REAL, subarray_type, ier)
diff --git a/src/generate_databases/parse_sep.c b/src/generate_databases/parse_sep.c
index fc7de14..a789c70 100644
--- a/src/generate_databases/parse_sep.c
+++ b/src/generate_databases/parse_sep.c
@@ -28,7 +28,7 @@
/*****************************************************************************/
void parse_sep_header(char *header_name,
int *n1, int *n2, int *n3,
- int *o1, int *o2, int *o3,
+ float *o1, float *o2, float *o3,
float *d1, float *d2, float *d3,
char *in) {
char *line = NULL;
@@ -36,7 +36,7 @@ void parse_sep_header(char *header_name,
FILE *fd = fopen(header_name, "r");
if (fd == NULL)
- EXIT_ON_ERR("fopen");
+ EXIT_ON_ERR(header_name);
/* Build regex mathing "name=value" with optional spaces*/
regex_t re;
@@ -68,11 +68,11 @@ void parse_sep_header(char *header_name,
} else if (!strcmp(name, "n3")) {
*n3 = atoi(value);
} else if (!strcmp(name, "o1")) {
- *o1 = atoi(value);
+ *o1 = atof(value);
} else if (!strcmp(name, "o2")) {
- *o2 = atoi(value);
+ *o2 = atof(value);
} else if (!strcmp(name, "o3")) {
- *o3 = atoi(value);
+ *o3 = atof(value);
} else if (!strcmp(name, "d1")) {
*d1 = atof(value);
} else if (!strcmp(name, "d2")) {
More information about the CIG-COMMITS
mailing list