From 65cb2a4cd938b6592e201fdb44f6b0558f70fce9 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Tue, 7 Mar 2023 15:17:53 +0300 Subject: [PATCH] updates user output for tomo models --- src/meshfem2D/meshfem2D.F90 | 317 +----------------- src/shared/parallel.F90 | 48 ++- src/shared/read_parameter_file.F90 | 14 +- .../define_external_model_from_tomo_file.f90 | 81 ++++- src/specfem2D/initialize_simulation.F90 | 2 +- src/specfem2D/read_mesh_databases.F90 | 2 +- src/specfem2D/setup_mesh.F90 | 6 +- src/specfem2D/specfem2D.F90 | 29 +- utils/createTomographyFile.py | 2 +- 9 files changed, 168 insertions(+), 333 deletions(-) diff --git a/src/meshfem2D/meshfem2D.F90 b/src/meshfem2D/meshfem2D.F90 index 079ad22e2..647205ed0 100644 --- a/src/meshfem2D/meshfem2D.F90 +++ b/src/meshfem2D/meshfem2D.F90 @@ -36,320 +36,10 @@ ! Basic mesh generator for SPECFEM2D ! !======================================================================== - -! If you use this code for your own research, please cite at least one article -! written by the developers of the package, for instance: -! -! @ARTICLE{TrKoLi08, -! author = {Jeroen Tromp and Dimitri Komatitsch and Qinya Liu}, -! title = {Spectral-Element and Adjoint Methods in Seismology}, -! journal = {communications in Computational Physics}, -! year = {2008}, -! volume = {3}, -! pages = {1-32}, -! number = {1}} -! -! @ARTICLE{PeKoLuMaLeCaLeMaLiBlNiBaTr11, -! author = {Daniel Peter and Dimitri Komatitsch and Yang Luo and Roland Martin -! and Nicolas {Le Goff} and Emanuele Casarotti and Pieyre {Le Loher} -! and Federica Magnoni and Qinya Liu and C\'eline Blitz and Tarje Nissen-Meyer -! and Piero Basini and Jeroen Tromp}, -! title = {Forward and adjoint simulations of seismic wave propagation on fully -! unstructured hexahedral meshes}, -! journal={Geophys. J. Int.}, -! year = {2011}, -! volume = {186}, -! pages = {721-739}, -! number = {2}, -! doi = {10.1111/j.1365-246X.2011.05044.x}} -! -! or -! -! @ARTICLE{VaCaSaKoVi99, -! author = {R. Vai and J. M. Castillo-Covarrubias and F. J. S\'anchez-Sesma and -! D. Komatitsch and J. P. Vilotte}, -! title = {Elastic wave propagation in an irregularly layered medium}, -! journal = {Soil Dynamics and Earthquake Engineering}, -! year = {1999}, -! volume = {18}, -! pages = {11-18}, -! number = {1}, -! doi = {10.1016/S0267-7261(98)00027-X}} -! -! @ARTICLE{LeChKoHuTr09, -! author = {Shiann Jong Lee and Yu Chang Chan and Dimitri Komatitsch and Bor -! Shouh Huang and Jeroen Tromp}, -! title = {Effects of realistic surface topography on seismic ground motion -! in the {Y}angminshan region of {T}aiwan based upon the spectral-element -! method and {LiDAR DTM}}, -! journal = {Bull. Seismol. Soc. Am.}, -! year = {2009}, -! volume = {99}, -! pages = {681-693}, -! number = {2A}, -! doi = {10.1785/0120080264}} -! -! @ARTICLE{LeChLiKoHuTr08, -! author = {Shiann Jong Lee and How Wei Chen and Qinya Liu and Dimitri Komatitsch -! and Bor Shouh Huang and Jeroen Tromp}, -! title = {Three-Dimensional Simulations of Seismic Wave Propagation in the -! {T}aipei Basin with Realistic Topography Based upon the Spectral-Element Method}, -! journal = {Bull. Seismol. Soc. Am.}, -! year = {2008}, -! volume = {98}, -! pages = {253-264}, -! number = {1}, -! doi = {10.1785/0120070033}} -! -! @ARTICLE{LeKoHuTr09, -! author = {S. J. Lee and Dimitri Komatitsch and B. S. Huang and J. Tromp}, -! title = {Effects of topography on seismic wave propagation: An example from -! northern {T}aiwan}, -! journal = {Bull. Seismol. Soc. Am.}, -! year = {2009}, -! volume = {99}, -! pages = {314-325}, -! number = {1}, -! doi = {10.1785/0120080020}} -! -! @ARTICLE{KoErGoMi10, -! author = {Dimitri Komatitsch and Gordon Erlebacher and Dominik G\"oddeke and -! David Mich\'ea}, -! title = {High-order finite-element seismic wave propagation modeling with -! {MPI} on a large {GPU} cluster}, -! journal = {J. Comput. Phys.}, -! year = {2010}, -! volume = {229}, -! pages = {7692-7714}, -! number = {20}, -! doi = {10.1016/j.jcp.2010.06.024}} -! -! @ARTICLE{KoGoErMi10, -! author = {Dimitri Komatitsch and Dominik G\"oddeke and Gordon Erlebacher and -! David Mich\'ea}, -! title = {Modeling the propagation of elastic waves using spectral elements -! on a cluster of 192 {GPU}s}, -! journal = {Computer Science Research and Development}, -! year = {2010}, -! volume = {25}, -! pages = {75-82}, -! number = {1-2}, -! doi = {10.1007/s00450-010-0109-1}} -! -! @ARTICLE{KoMiEr09, -! author = {Dimitri Komatitsch and David Mich\'ea and Gordon Erlebacher}, -! title = {Porting a high-order finite-element earthquake modeling application -! to {NVIDIA} graphics cards using {CUDA}}, -! journal = {Journal of Parallel and Distributed Computing}, -! year = {2009}, -! volume = {69}, -! pages = {451-460}, -! number = {5}, -! doi = {10.1016/j.jpdc.2009.01.006}} -! -! @ARTICLE{LiPoKoTr04, -! author = {Qinya Liu and Jascha Polet and Dimitri Komatitsch and Jeroen Tromp}, -! title = {Spectral-element moment tensor inversions for earthquakes in {S}outhern {C}alifornia}, -! journal={Bull. Seismol. Soc. Am.}, -! year = {2004}, -! volume = {94}, -! pages = {1748-1761}, -! number = {5}, -! doi = {10.1785/012004038}} -! -! @INCOLLECTION{ChKoViCaVaFe07, -! author = {Emmanuel Chaljub and Dimitri Komatitsch and Jean-Pierre Vilotte and -! Yann Capdeville and Bernard Valette and Gaetano Festa}, -! title = {Spectral Element Analysis in Seismology}, -! booktitle = {Advances in Wave Propagation in Heterogeneous Media}, -! publisher = {Elsevier - Academic Press}, -! year = {2007}, -! editor = {Ru-Shan Wu and Val\'erie Maupin}, -! volume = {48}, -! series = {Advances in Geophysics}, -! pages = {365-419}} -! -! @ARTICLE{KoVi98, -! author={D. Komatitsch and J. P. Vilotte}, -! title={The spectral-element method: an efficient tool to simulate the seismic response of 2{D} and 3{D} geological structures}, -! journal={Bull. Seismol. Soc. Am.}, -! year=1998, -! volume=88, -! number= 2, -! pages={368-392}} -! -! @ARTICLE{KoTr99, -! author={D. Komatitsch and J. Tromp}, -! year=1999, -! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation}, -! journal={Geophys. J. Int.}, -! volume=139, -! number=3, -! pages={806-822}, -! doi={10.1046/j.1365-246x.1999.00967.x}} -! -! @ARTICLE{KoLiTrSuStSh04, -! author={Dimitri Komatitsch and Qinya Liu and Jeroen Tromp and Peter S\"{u}ss -! and Christiane Stidham and John H. Shaw}, -! year=2004, -! title={Simulations of Ground Motion in the {L}os {A}ngeles {B}asin -! based upon the Spectral-Element Method}, -! journal={Bull. Seism. Soc. Am.}, -! volume=94, -! number= 1, -! pages={187-206}} -! -! @ARTICLE{MoTr08, -! author={C. Morency and J. Tromp}, -! title={Spectral-element simulations of wave propagation in poroelastic media}, -! journal={Geophys. J. Int.}, -! year=2008, -! volume=175, -! pages={301-345}} -! -! and/or other articles from https://specfem.org/komatitsch.free.fr/publications.html ! -! If you use the kernel capabilities of the code, please cite at least one article -! written by the developers of the package, for instance: +! Please find in the header of specfem2D.F90 further code informations. ! -! @ARTICLE{TrKoLi08, -! author = {Jeroen Tromp and Dimitri Komatitsch and Qinya Liu}, -! title = {Spectral-Element and Adjoint Methods in Seismology}, -! journal = {communications in Computational Physics}, -! year = {2008}, -! volume = {3}, -! pages = {1-32}, -! number = {1}} -! -! @ARTICLE{PeKoLuMaLeCaLeMaLiBlNiBaTr11, -! author = {Daniel Peter and Dimitri Komatitsch and Yang Luo and Roland Martin -! and Nicolas {Le Goff} and Emanuele Casarotti and Pieyre {Le Loher} -! and Federica Magnoni and Qinya Liu and C\'eline Blitz and Tarje Nissen-Meyer -! and Piero Basini and Jeroen Tromp}, -! title = {Forward and adjoint simulations of seismic wave propagation on fully -! unstructured hexahedral meshes}, -! journal={Geophys. J. Int.}, -! year = {2011}, -! volume = {186}, -! pages = {721-739}, -! number = {2}, -! doi = {10.1111/j.1365-246X.2011.05044.x}} -! -! @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}, -! title={Finite-frequency kernels for wave propagation in porous media based upon adjoint methods}, -! year=2009, -! journal={Geophys. J. Int.}, -! doi={10.1111/j.1365-246X.2009.04332}} -! -! If you use the METIS / SCOTCH / CUBIT non-structured capabilities, please also cite: -! -! @ARTICLE{MaKoBlLe08, -! author = {R. Martin and D. Komatitsch and C. Blitz and N. {Le Goff}}, -! title = {Simulation of seismic wave propagation in an asteroid based upon -! an unstructured {MPI} spectral-element method: blocking and non-blocking -! communication strategies}, -! journal = {Lecture Notes in Computer Science}, -! year = {2008}, -! volume = {5336}, -! pages = {350-363}} -! -! -! version 8.0, Etienne Bachmann, Alexis Bottero, Bryant Chow, Paul Cristini, Rene Gassmoeller, Michael Gineste, -! Felix Halpaap, Dimitri Komatitsch, Matthieu Lefebvre, Qiancheng Liu, Qinya Liu, Zhaolun Liu, -! David Luet, Ryan Modrak, Christina Morency, Daniel Peter, Eric Rosenkrantz, Herurisa Rusmanugroho, -! Elliott Sales de Andrade, Eduardo Valero Cano, Zhinan Xie, Zhendong Zhang, December 2022: -! - various code improvements -! - GPU support -! - axisymmetric 2.5D simulations -! -! version 7.0, Dimitri Komatitsch, Zhinan Xie, Paul Cristini, Roland Martin and Rene Matzen, July 2012: -! - added support for Convolution PML absorbing layers -! - added higher-order time schemes (4th order Runge-Kutta and LDDRK4-6) -! - many small or moderate bug fixes -! -! version 6.2, many developers, April 2011: -! - restructured package source code into separate src/ directories -! - added configure & Makefile scripts and a PDF manual in doc/ -! - added user examples in EXAMPLES/ -! - added a USER_T0 parameter to fix the onset time in simulation -! -! version 6.1, Christina Morency and Pieyre Le Loher, March 2010: -! - added SH (membrane) waves calculation for elastic media -! - added support for external fully anisotropic media -! - fixed some bugs in acoustic kernels -! -! version 6.0, Christina Morency and Yang Luo, August 2009: -! - support for poroelastic media -! - adjoint method for acoustic/elastic/poroelastic -! -! version 5.2, Dimitri Komatitsch, Nicolas Le Goff and Roland Martin, February 2008: -! - support for CUBIT and GiD meshes -! - MPI implementation of the code based on domain decomposition -! with METIS or SCOTCH -! - general fluid/solid implementation with any number, shape and orientation of -! matching edges -! - fluid potential of density * displacement instead of displacement -! - absorbing edges with any normal vector -! - general numbering of absorbing and acoustic free surface edges -! - cleaned implementation of attenuation as in Carcione (1993) -! - merged loops in the solver for efficiency -! - simplified input of external model -! - added CPU time information -! - translated many comments from French to English -! -! version 5.1, Dimitri Komatitsch, January 2005: -! - more general mesher with any number of curved layers -! - Dirac and Gaussian time sources and corresponding convolution routine -! - option for acoustic medium instead of elastic -! - receivers at any location, not only grid points -! - moment-tensor source at any location, not only a grid point -! - color snapshots -! - more flexible DATA/Par_file with any number of comment lines -! - Xsu scripts for seismograms -! - subtract t0 from seismograms -! - seismograms and snapshots in pressure in addition to vector field -! -! version 5.0, Dimitri Komatitsch, May 2004: -! - got rid of useless routines, suppressed commons etc. -! - weak formulation based explicitly on stress tensor -! - implementation of full anisotropy -! - implementation of attenuation based on memory variables -! -! based on SPECFEM2D version 4.2, June 1998 -! (c) by Dimitri Komatitsch, Harvard University, USA -! and Jean-Pierre Vilotte, Institut de Physique du Globe de Paris, France -! -! itself based on SPECFEM2D version 1.0, 1995 -! (c) by Dimitri Komatitsch and Jean-Pierre Vilotte, -! Institut de Physique du Globe de Paris, France -! - -! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette, -! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential -! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002). -! This permits acoustic-elastic coupling based on a non-iterative time scheme. -! Displacement is then: u = grad(Chi) / rho -! Velocity is then: v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi) -! and pressure is: p = - Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi). -! The source in an acoustic element is a pressure source. -! First-order acoustic-acoustic discontinuities are also handled automatically -! because pressure is continuous at such an interface, therefore Chi_dot_dot -! is continuous, therefore Chi is also continuous, which is consistent with -! the spectral-element basis functions and with the assembling process. -! This is the reason why a simple displacement potential u = grad(Chi) would -! not work because it would be discontinuous at such an interface and would -! therefore not be consistent with the basis functions. +! ************** PROGRAM STARTS HERE ************** program meshfem2D @@ -357,7 +47,6 @@ program meshfem2D use shared_parameters use part_unstruct_par -! use source_file_par use compute_elements_load_par implicit none @@ -412,7 +101,7 @@ program meshfem2D ! *** ! reads in parameters in DATA/Par_file BROADCAST_AFTER_READ = .false. - call read_parameter_file(1,BROADCAST_AFTER_READ) + call read_parameter_file(.true.,BROADCAST_AFTER_READ) ! reads in additional files for mesh elements if (read_external_mesh) then diff --git a/src/shared/parallel.F90 b/src/shared/parallel.F90 index cd357be64..d9d497b99 100644 --- a/src/shared/parallel.F90 +++ b/src/shared/parallel.F90 @@ -691,7 +691,7 @@ subroutine min_all_i(sendbuf, recvbuf) implicit none - integer:: sendbuf, recvbuf + integer :: sendbuf, recvbuf #ifdef WITH_MPI ! local parameters @@ -727,6 +727,52 @@ subroutine max_all_i(sendbuf, recvbuf) end subroutine max_all_i +! +!------------------------------------------------------------------------------------------------- +! + + subroutine min_all_dp(sendbuf, recvbuf) + + use my_mpi + + implicit none + + double precision :: sendbuf, recvbuf + +#ifdef WITH_MPI + ! local parameters + integer ier + + call MPI_REDUCE(sendbuf,recvbuf,1,MPI_DOUBLE_PRECISION,MPI_MIN,0,my_local_mpi_comm_world,ier) +#else + recvbuf = sendbuf +#endif + + end subroutine min_all_dp + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine max_all_dp(sendbuf, recvbuf) + + use my_mpi + + implicit none + + double precision :: sendbuf, recvbuf + +#ifdef WITH_MPI + ! local parameters + integer :: ier + + call MPI_REDUCE(sendbuf,recvbuf,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,my_local_mpi_comm_world,ier) +#else + recvbuf = sendbuf +#endif + + end subroutine max_all_dp + ! !------------------------------------------------------------------------------------------------- ! diff --git a/src/shared/read_parameter_file.F90 b/src/shared/read_parameter_file.F90 index 5e36f9d7b..389c87b20 100644 --- a/src/shared/read_parameter_file.F90 +++ b/src/shared/read_parameter_file.F90 @@ -31,7 +31,7 @@ ! !======================================================================== - subroutine read_parameter_file(imesher,BROADCAST_AFTER_READ) + subroutine read_parameter_file(is_mesher,BROADCAST_AFTER_READ) ! reads in DATA/Par_file @@ -40,7 +40,7 @@ subroutine read_parameter_file(imesher,BROADCAST_AFTER_READ) implicit none - integer, intent(in) :: imesher + logical, intent(in) :: is_mesher logical, intent(in) :: BROADCAST_AFTER_READ ! initializes @@ -63,12 +63,14 @@ subroutine read_parameter_file(imesher,BROADCAST_AFTER_READ) ! user output write(IMAIN,*) 'Title of the simulation: ',trim(title) write(IMAIN,*) - if (AXISYM) write(IMAIN,*) 'Axisymmetric simulation' - write(IMAIN,*) + if (AXISYM) then + write(IMAIN,*) 'Axisymmetric simulation' + write(IMAIN,*) + endif call flush_IMAIN() ! reads receiver lines - if (imesher == 1) then + if (is_mesher) then ! user output write(IMAIN,*) 'Receiver lines:' write(IMAIN,*) ' Nb of line sets = ',nreceiversets @@ -78,7 +80,7 @@ subroutine read_parameter_file(imesher,BROADCAST_AFTER_READ) call read_parameter_file_receiversets() ! only mesher needs to reads this - if (imesher == 1) then + if (is_mesher) then ! reads material definitions call read_material_table() diff --git a/src/specfem2D/define_external_model_from_tomo_file.f90 b/src/specfem2D/define_external_model_from_tomo_file.f90 index 947808d77..870d73651 100644 --- a/src/specfem2D/define_external_model_from_tomo_file.f90 +++ b/src/specfem2D/define_external_model_from_tomo_file.f90 @@ -365,6 +365,8 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & ! local parameters integer :: i,j,ispec,iglob double precision :: xmesh,zmesh + ! stats + integer :: npoint_tomo,npoint_internal ! user output if (myrank == 0) then @@ -392,6 +394,8 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & c22ext(:,:,:) = 0.d0 ! for AXISYM only ! loop on all the elements of the mesh, and inside each element loop on all the GLL points + npoint_tomo = 0 + npoint_internal = 0 do ispec = 1,nspec do j = 1,NGLLZ do i = 1,NGLLX @@ -425,6 +429,10 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & QKappa_attenuationcoef(kmato(ispec)) = QKappa_attenuationext(i,j,ispec) Qmu_attenuationcoef(kmato(ispec)) = Qmu_attenuationext(i,j,ispec) + + ! counting gll points with tomography model values assigned + npoint_tomo = npoint_tomo + 1 + else ! Internal model rhoext(i,j,ispec) = density(1,kmato(ispec)) @@ -444,11 +452,32 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & c23ext(i,j,ispec) = anisotropycoef(8,kmato(ispec)) c25ext(i,j,ispec) = anisotropycoef(9,kmato(ispec)) c22ext(i,j,ispec) = anisotropycoef(10,kmato(ispec)) ! for AXISYM + + ! counting gll points with internal model values assigned + npoint_tomo = npoint_tomo + 1 endif enddo enddo enddo + call synchronize_all() + + ! collect totals + iglob = npoint_tomo + call max_all_i(iglob,npoint_tomo) + iglob = npoint_internal + call max_all_i(iglob,npoint_internal) + + ! user output + if (myrank == 0) then + if (npoint_tomo > 0) & + write(IMAIN,*) ' number of GLL points with tomographic model values: ',npoint_tomo + if (npoint_internal > 0) & + write(IMAIN,*) ' number of GLL points with internal model values : ',npoint_internal + write(IMAIN,*) + call flush_IMAIN() + endif + end subroutine define_external_model_from_tomo_file ! @@ -509,7 +538,7 @@ subroutine read_tomo_file() use specfem_par, only: myrank,TOMOGRAPHY_FILE use model_tomography_par - use constants, only: IIN,IMAIN + use constants, only: IIN,IMAIN,HUGEVAL implicit none @@ -518,6 +547,7 @@ subroutine read_tomo_file() character(len=150) :: string_read double precision, dimension(:), allocatable :: x_tomography,z_tomography,vp_tomography,vs_tomography,rho_tomography + double precision :: read_rho_min,read_rho_max,read_vp_min,read_vp_max,read_vs_min,read_vs_max,tmp_dp ! opens file for reading open(unit=IIN,file=trim(TOMOGRAPHY_FILE),status='old',action='read',iostat=ier) @@ -562,11 +592,26 @@ subroutine read_tomo_file() ! allocate models parameter records allocate(x_tomography(nrecord),z_tomography(nrecord),vp_tomography(nrecord),vs_tomography(nrecord), & rho_tomography(nrecord),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate tomo arrays') + x_tomography(:) = 0.d0; z_tomography(:) = 0.d0 + vp_tomography(:) = 0.d0; vs_tomography(:) = 0.d0; rho_tomography(:) = 0.d0 + allocate(x_tomo(NX),z_tomo(NZ),vp_tomo(NX,NZ),vs_tomo(NX,NZ),rho_tomo(NX,NZ),stat=ier) if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate tomo arrays') + x_tomo(:) = 0.d0; z_tomo(:) = 0.d0 + vp_tomo(:,:) = 0.d0; vs_tomo(:,:) = 0.d0; rho_tomo(:,:) = 0.d0 ! Checks the number of records for points definition while storing them irecord = 0 + + ! stats + read_vp_min = + HUGEVAL + read_vp_max = - HUGEVAL + read_vs_min = + HUGEVAL + read_vs_max = - HUGEVAL + read_rho_min = + HUGEVAL + read_rho_max = - HUGEVAL + do while (irecord < nrecord) call tomo_read_next_line(IIN,string_read) read(string_read,*) x_tomography(irecord+1),z_tomography(irecord+1), & @@ -574,6 +619,15 @@ subroutine read_tomo_file() if (irecord < NX) x_tomo(irecord+1) = x_tomography(irecord+1) + ! stats + if (vp_tomography(irecord+1) > read_vp_max) read_vp_max = vp_tomography(irecord+1) + if (vp_tomography(irecord+1) < read_vp_min) read_vp_min = vp_tomography(irecord+1) + if (vs_tomography(irecord+1) > read_vs_max) read_vs_max = vs_tomography(irecord+1) + if (vs_tomography(irecord+1) < read_vs_min) read_vs_min = vs_tomography(irecord+1) + if (rho_tomography(irecord+1) > read_rho_max) read_rho_max = rho_tomography(irecord+1) + if (rho_tomography(irecord+1) < read_rho_min) read_rho_min = rho_tomography(irecord+1) + + ! counter irecord = irecord + 1 enddo @@ -587,16 +641,39 @@ subroutine read_tomo_file() enddo enddo + call synchronize_all() + if (irecord /= nrecord .and. myrank == 0) then print *, 'Error: ',trim(TOMOGRAPHY_FILE),' has invalid number of records' - print *, ' number of grid points specified (= NX*NZ):',nrecord + print *, ' number of grid points specified (= NX*NZ) :',nrecord print *, ' number of file lines for grid points :',irecord call stop_the_code('Error in tomography data file for the grid points definition') endif + ! collects stats + tmp_dp = read_vp_min + call min_all_dp(tmp_dp,read_vp_min) + tmp_dp = read_vp_max + call max_all_dp(tmp_dp,read_vp_max) + + tmp_dp = read_vs_min + call min_all_dp(tmp_dp,read_vs_min) + tmp_dp = read_vs_max + call max_all_dp(tmp_dp,read_vs_max) + + tmp_dp = read_rho_min + call min_all_dp(tmp_dp,read_rho_min) + tmp_dp = read_rho_max + call max_all_dp(tmp_dp,read_rho_max) + + ! user output if (myrank == 0) then write(IMAIN,*) ' Number of grid points = NX*NZ:',nrecord write(IMAIN,*) + write(IMAIN,*) ' read model vp : min/max = ',sngl(read_vp_min),' / ',sngl(read_vp_max) + write(IMAIN,*) ' vs : min/max = ',sngl(read_vs_min),' / ',sngl(read_vs_max) + write(IMAIN,*) ' rho: min/max = ',sngl(read_rho_min),' / ',sngl(read_rho_max) + write(IMAIN,*) call flush_IMAIN() endif diff --git a/src/specfem2D/initialize_simulation.F90 b/src/specfem2D/initialize_simulation.F90 index b1d6e1abb..ba3ab3d57 100644 --- a/src/specfem2D/initialize_simulation.F90 +++ b/src/specfem2D/initialize_simulation.F90 @@ -128,7 +128,7 @@ subroutine initialize_simulation() ! read the parameter file BROADCAST_AFTER_READ = .true. - call read_parameter_file(0,BROADCAST_AFTER_READ) + call read_parameter_file(.false.,BROADCAST_AFTER_READ) ! reads in source descriptions ! note: we will need a source frequency for outputting poroelastic velocities when reading the mesh databases diff --git a/src/specfem2D/read_mesh_databases.F90 b/src/specfem2D/read_mesh_databases.F90 index 967516160..bde5a6e2d 100644 --- a/src/specfem2D/read_mesh_databases.F90 +++ b/src/specfem2D/read_mesh_databases.F90 @@ -1305,6 +1305,7 @@ subroutine read_mesh_databases_free_surf() if (myrank == 0) then write(IMAIN,*) 'reading free surface information...' write(IMAIN,*) ' number of acoustic free surface boundary elements = ',nelem_acoustic_surface + write(IMAIN,*) call flush_IMAIN() endif @@ -1345,7 +1346,6 @@ subroutine read_mesh_databases_free_surf() if (nelem_acoustic_surface_all > 0) then ! user output if (myrank == 0) then - write(IMAIN,*) write(IMAIN,*) ' Total number of acoustic free surface elements: ',nelem_acoustic_surface_all write(IMAIN,*) call flush_IMAIN() diff --git a/src/specfem2D/setup_mesh.F90 b/src/specfem2D/setup_mesh.F90 index 0ad4d0b72..8a0e79765 100644 --- a/src/specfem2D/setup_mesh.F90 +++ b/src/specfem2D/setup_mesh.F90 @@ -605,7 +605,11 @@ subroutine setup_mesh_material_properties() if (assign_external_model) then ! user output if (myrank == 0) then - write(IMAIN,*) ' assigning an external velocity and density model' + if (trim(MODEL) == 'tomo') then + write(IMAIN,*) ' assigning an external tomography velocity and density model' + else + write(IMAIN,*) ' assigning an external velocity and density model' + endif call flush_IMAIN() endif diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90 index c8659d2d1..afe3a509d 100644 --- a/src/specfem2D/specfem2D.F90 +++ b/src/specfem2D/specfem2D.F90 @@ -38,7 +38,7 @@ ! for the anelastic anisotropic or poroelastic wave equation. ! !==================================================================================== - +! ! If you use this code for your own research, please cite at least one article ! written by the developers of the package, for instance: ! @@ -266,6 +266,13 @@ ! volume = {5336}, ! pages = {350-363}} ! +! +! To report bugs or suggest improvements to the code, please use our online +! Issues tracking system at https://github.com/SPECFEM/specfem2d/ . +! +! Evolution of the code: +! --------------------- +! ! version 8.0, Etienne Bachmann, Alexis Bottero, Bryant Chow, Paul Cristini, Rene Gassmoeller, Michael Gineste, ! Felix Halpaap, Dimitri Komatitsch, Matthieu Lefebvre, Qiancheng Liu, Qinya Liu, Zhaolun Liu, ! David Luet, Ryan Modrak, Christina Morency, Daniel Peter, Eric Rosenkrantz, Herurisa Rusmanugroho, @@ -335,15 +342,23 @@ ! (c) by Dimitri Komatitsch and Jean-Pierre Vilotte, ! Institut de Physique du Globe de Paris, France ! - -! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette, +! +! In case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette, ! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential ! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002). +! ! This permits acoustic-elastic coupling based on a non-iterative time scheme. -! Displacement is then: u = grad(Chi) / rho -! Velocity is then: v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi) -! and pressure is: p = - Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi). +! Displacement is then: +! u = grad(Chi) / rho +! Velocity is then: +! v = grad(Chi_dot) / rho +! (Chi_dot being the time derivative of Chi) +! and pressure is: +! p = - Chi_dot_dot +! (Chi_dot_dot being the time second derivative of Chi). +! ! The source in an acoustic element is a pressure source. +! ! First-order acoustic-acoustic discontinuities are also handled automatically ! because pressure is continuous at such an interface, therefore Chi_dot_dot ! is continuous, therefore Chi is also continuous, which is consistent with @@ -351,6 +366,8 @@ ! This is the reason why a simple displacement potential u = grad(Chi) would ! not work because it would be discontinuous at such an interface and would ! therefore not be consistent with the basis functions. +! +! ************** PROGRAM STARTS HERE ************** program specfem2D diff --git a/utils/createTomographyFile.py b/utils/createTomographyFile.py index 034b9a7b3..d5a897482 100755 --- a/utils/createTomographyFile.py +++ b/utils/createTomographyFile.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/env python # -*- coding: utf-8 -*- """ Example of external velocity model