Commit 1df104fb authored by laura.gassner's avatar laura.gassner

visco-acoustic modelling is now available - not thoroughly tested yet.

parent 67c207a1
......@@ -248,7 +248,7 @@ In these files, each material parameter value must be saved as 32 bit (4 byte) n
\textit{ximage n1=$<$NY$>$ $<$ model/test.vp} .
\newline
It is also possible to read Qp, and Qs grid files to allow for spatial variable attenuation. For this you must uncomment a few lines in \texttt{src/readmod.c} and generate the corresponding binary files.
It is also possible to read Qp, and Qs grid files to allow for spatial variable attenuation. If no models are provided, the value TAU is used to generate Q-Models (see section \ref{q_approx}).
If READMOD=0 the model is generated ''on the fly'' by DENISE, i.e. it is generated internally before the time loop starts. See \texttt{genmod/1D\_linear\_gradient\_el.c} for an example function that generates a simple model with a linear vertical gradient ''on the fly''. If READMOD=0 this function is called in \texttt{src/denise.c} and therefore must be specified in \texttt{src/Makefile} (at the top of \texttt{src/Makefile}, see section \ref{compexec}). If you change this file, for example to change the model structure, you need to re-compile DENISE by changing to the src directory and ''make denise''.
......@@ -342,8 +342,10 @@ Default values are:
If SEISMO$>$0 seismograms recorded at the receiver positions are written to the corresponding output files. The sampling rate of the seismograms is NDT*DT seconds. In case of a small
time step interval and a high number of time steps, it might be useful to choose a high NDT in order to avoid a unnecessary detailed sampling of the seismograms and consequently large files of seismogram data. Possible output formats of the seismograms are SU, ASCII and BINARY. It is recommended to use SU format for saving the seismograms. The main advantage of this format is that the time step interval (NDT*DT) and the acquisition geometry (shot and receiver locations) are stored in the corresponding SU header words. Also additional header words like offset are set by DENISE. This format thus facilitates a further visualization and processing of the synthetic seismograms. Note, however, that SU cannot handle sampling rates smaller than 1.0e-6 seconds and the number of samples is limited to about 32.000. In such cases, you should increase the sampling rate by increasing NDT. If this is impossible (for example because the Nyquist criterion is violated) you must choose a different output format (ASCII or binary).
\newpage
\section{Q-approximation}
\label{q_approx}
{\color{blue}{\begin{verbatim}
"Q-approximation",
"L" : "0",
......@@ -724,6 +726,7 @@ With VP\_VS\_RATIO = 1.5 (e.g.) it is possible to ensure a minimum Vp/Vs ratio d
{\color{blue}{\begin{verbatim}
"Time windowing and damping" : "comment",
"TIMEWIN" : "0",
"TW_IND" : "0",
"PICKS_FILE" : "./picked_times/picks"
"TWLENGTH_PLUS" : "0.01",
"TWLENGTH_MINUS" : "0.01",
......@@ -736,7 +739,7 @@ Default values are:
\end{verbatim}}}
To apply time windowing in a time series the paramter TIMEWIN must set to 1. A automatic picker routine is not integrated at the moment. The point in time (picked time) for each source must be specified in seperate files. The folder and file name can be set with the parameter PICKS\_FILE. The files must be named like this PICKS\_FILE\_<sourcenumber>.dat. So the number of sources in (SRCREC) must be equal to the number of files. Each file must contain the picked times for every receiver.\
The parameters TWLENGTH\_PLUS and TWLENGTH\_MINUS specify the length of the time window after (PLUS) and before (MINUS) the picked time. The unit is seconds. The damping factor GAMMA must be set individually.
The parameters TWLENGTH\_PLUS and TWLENGTH\_MINUS specify the length of the time window after (PLUS) and before (MINUS) the picked time. The unit is seconds. The damping factor GAMMA must be set individually. When usig TW\_IND = 1 three columns are expected in PICK\_FILE for each trace with the picked time, time window length plus and time window length minus. In this case TWLENGTH\_PLUS and TWLENGTH\_MINUS are ignored.
\section{Source wavelet inversion}
......@@ -749,6 +752,7 @@ it is necessary to design a filter which minimizes the misfit to the field recor
"PARA" : "fdlsq:tshift=0.0",
"N_STF" : "10",
"N_STF_START" : "1",
"TAPER_STF" : "0",
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill.dat",
\end{verbatim}}}
......@@ -769,7 +773,9 @@ Examples for the parameter PARA are:
\textit{fdlsq:exp=1.4:pow2}
\end{itemize}
N\_STF is the increment between the iteration steps. N\_STF\_START defines at which iterationstep the inversion for STF should start. This parameter has to be set at least to 1 NOT(!) 0. With TRKILL\_STF = 1 it is possible to apply a trace killing for the estimation of the source wavelet correction filter.
N\_STF is the increment between the iteration steps. N\_STF\_START defines at which iterationstep the inversion for STF should start. This parameter has to be set at least to 1 NOT(!) 0.
When using TAPER\_STF = 1, the source signal is tapered. See \texttt{src/taper.c} for the taper definition.
With TRKILL\_STF = 1 it is possible to apply a trace killing for the estimation of the source wavelet correction filter.
\newline
Please note: If you additionally switch on frequency filtering during the inversion (TIME\_FILT=1 or TIME\_FILT=2), the parameters N\_STF and N\_STF\_START will be ignored. But the optimal source time function will be inverted for the first iteration and after every change of the frequency range.
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2005 <name of author>
*
* This file is part of DENISE.
*
* DENISE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* DENISE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*
* Model homogeneous half space
* last update 11.04.02, T. Bohlen
*/
#include "fd.h"
void model_viscac(float ** rho, float ** pi, float ** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NXG, NYG, POS[3], L, MYID;
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
extern float DH, *FL, TAU, DT;
/* local variables */
float vp, rhov, grad1, grad3, y, tp, piv, *pts;
int i, j, ii, jj, l;
char modfile[STRING_SIZE];
/* parameters for layer 1 velocity in m/s, depth h in meter */
const float vp1=600.0, rho1=1800.0, h=10.0;
/* parameters for layer 2 due to calculation of grad1, grad2 and grad3*/
const float vp2=2500.0, rho2=2000.0;
/*-----------------------------------------------------------------------*/
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
y=h/DH;
if(y==NYG) err(" \n y is equal NYG !! see src/model_grad.c \n ");
grad1=(vp2-vp1)/y;
grad3=(rho2-rho1)/y;
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
if(j<=y){
vp=vp1+(j*grad1);
rhov=rho1+(j*grad3);
}
else{
vp=vp2;
rhov=rho2;
}
tp=TAU;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
rho[jj][ii]=rhov;
pi[jj][ii]=vp;
taup[jj][ii]=tp;
}
}
}
sprintf(modfile,"%s_rho_it_0.bin",INV_MODELFILE);
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
sprintf(modfile,"%s_vp_it_0.bin",INV_MODELFILE);
writemod(modfile,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
free_vector(pts,1,L);
}
......@@ -7,8 +7,9 @@
#MODEL = hh.c
MODEL = ../genmod/1D_linear_gradient_visc.c
MODEL_AC = ../genmod/1D_linear_gradient_ac.c
MODEL_EL = ../genmod/1D_linear_gradient_el.c
MODEL_AC = ../genmod/1D_linear_gradient_ac.c
MODEL_VAC = ../genmod/1D_linear_gradient_viscac.c
EXEC= ../bin
# Compiler (LAM: CC=hcc, CRAY T3E: CC=cc)
......@@ -82,7 +83,9 @@ DENISE= \
update_v_PML.c \
update_v_acoustic_PML.c \
prepare_update_s.c \
prepare_update_p.c \
update_p_PML.c \
update_p_visc_PML.c \
update_s_elastic_ssg.c \
update_s_elastic_PML.c \
update_s_visc_PML.c \
......@@ -101,11 +104,13 @@ DENISE= \
LBFGS1.c \
smooth.c \
$(MODEL) \
$(MODEL_AC) \
$(MODEL_EL) \
$(MODEL_AC) \
$(MODEL_VAC) \
matcopy.c \
matcopy_elastic.c \
matcopy_acoustic.c \
matcopy_viscac.c \
mergemod.c \
max_grad.c \
note.c \
......@@ -123,6 +128,7 @@ DENISE= \
readmod.c \
readmod_elastic.c \
readmod_acoustic.c \
readmod_viscac.c \
receiver.c \
rd_sour.c \
saveseis.c \
......@@ -151,6 +157,7 @@ DENISE= \
zero_fdveps.c \
zero_fdveps_ac.c \
zero_fdveps_visc.c \
zero_fdveps_viscac.c \
calc_envelope.c \
calc_hilbert.c \
filter_frequencies.c
......
This diff is collapsed.
......@@ -148,26 +148,20 @@ void LBFGS1(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos,
double LU_decomp(double **A, double *x, double *b,int n);
float minimum_m(float **mat, int nx, int ny);
float maximum_m(float **mat, int nx, int ny);
float maximum_m(float **mat, int nx, int ny);
void model(float ** rho, float ** pi, float ** u,
float ** taus, float ** taup, float * eta);
void model(float ** rho, float ** pi, float ** u, float ** taus, float ** taup, float * eta);
void model_elastic(float ** rho, float ** pi, float ** u);
void model_ani(float ** rho, float ** c11, float ** c15, float ** c13,
float ** c35, float ** c33, float ** c55,
float ** taus, float ** taup, float * eta);
void model_ani(float ** rho, float ** c11, float ** c15, float ** c13, float ** c35, float ** c33, float ** c55, float ** taus, float ** taup, float * eta);
void matcopy(float ** prho, float ** ppi, float ** pu, float ** ptaup,
float ** ptaus);
void matcopy(float ** prho, float ** ppi, float ** pu, float ** ptaup, float ** ptaus);
void matcopy_elastic(float ** prho, float ** ppi, float ** pu);
void matcopy_ani(float ** rho, float ** c11, float ** c15, float ** c13,
float ** c35, float ** c33, float ** c55, float ** taus,
float ** taup);
void matcopy_ani(float ** rho, float ** c11, float ** c15, float ** c13, float ** c35, float ** c33, float ** c55, float ** taus, float ** taup);
void max_grad(float ** waveconv, float ** waveconv_rho, float ** waveconv_u, float ** rho, float ** pi, float ** u);
......@@ -231,8 +225,7 @@ float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy);
void read_par_json(FILE *fp, char *fileinp);
void readmod(float ** rho, float ** pi, float ** u,
float ** taus, float ** taup, float * eta);
void readmod(float ** rho, float ** pi, float ** u, float ** taus, float ** taup, float * eta);
void readmod_elastic(float ** rho, float ** pi, float ** u);
......@@ -498,7 +491,6 @@ void zero(float *A, int u_max);
void normalize_data(float **data, int ntr, int ns);
/* functions for acoustic modelling */
void model_acoustic(float ** rho, float ** pi);
void readmod_acoustic(float ** rho, float ** pi);
void matcopy_acoustic(float ** prho, float ** ppi);
......@@ -506,22 +498,31 @@ void zero_fdveps_ac(int ny1, int ny2, int nx1, int nx2, float ** vx, float ** vy
float ** psi_sxx_x, float ** psi_sxy_x, float ** psi_vxx, float ** psi_vyx,
float ** psi_syy_y, float ** psi_sxy_y, float ** psi_vyy, float ** psi_vxy, float ** psi_vxxs);
void update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nt,
float ** vx, float ** vxp1, float ** vxm1, float ** vy, float ** vyp1, float ** vym1,
float ** sp, float **rip, float **rjp, float ** srcpos_loc, float ** signals, float ** signals1, int nsrc, float ** absorb_coeff,
float *hc, int infoout,int sw, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_sxx_x, float ** psi_syy_y);
void update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nt, float ** vx, float ** vxp1, float ** vxm1, float ** vy, float ** vyp1, float ** vym1,
float ** sp, float **rip, float **rjp, float ** srcpos_loc, float ** signals, float ** signals1, int nsrc, float ** absorb_coeff, float *hc, int infoout, int sw,
float * K_x_half, float * a_x_half, float * b_x_half, float * K_y_half, float * a_y_half, float * b_y_half, float ** psi_sxx_x, float ** psi_syy_y);
void update_p_PML(int nx1, int nx2, int ny1, int ny2,
float ** vx, float ** vy, float ** sp, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx);
void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sp, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx);
void surface_acoustic_PML(int ndepth, float ** sp);
void exchange_p(float ** sp, float ** bufferlef_to_rig, float ** bufferrig_to_lef,
float ** buffertop_to_bot, float ** bufferbot_to_top,
MPI_Request * req_send, MPI_Request * req_rec);
void exchange_p(float ** sp, float ** bufferlef_to_rig, float ** bufferrig_to_lef, float ** buffertop_to_bot, float ** bufferbot_to_top, MPI_Request * req_send, MPI_Request * req_rec);
/* functions for viscoacoustic modelling */
void model_viscac(float ** rho, float ** pi, float ** taup, float * eta);
void readmod_viscac(float ** rho, float ** pi, float ** taup, float * eta);
void matcopy_viscac(float ** prho, float ** ppi, float ** taup);
void prepare_update_p(float *etajm, float *peta, float **ppi, float **prho, float **ptaup, float **g, float *bjm, float *cjm, float ***e);
void zero_fdveps_viscac(int ny1, int ny2, int nx1, int nx2, float ** vx, float ** vy, float ** sp, float ** vxp1, float ** vyp1,
float ** psi_sxx_x, float ** psi_sxy_x, float ** psi_vxx, float ** psi_vyx, float ** psi_syy_y, float ** psi_sxy_y, float ** psi_vyy, float ** psi_vxy, float ** psi_vxxs, float ***pp);
void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sp, float ** pi, float **rho, float *hc, int infoout,
float ***p, float **g, float *bjm, float *cjm, float ***e,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx);
\ No newline at end of file
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* DENISE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* For the averaging of material properties each process requires values
* at the indices 0 and NX+1 etc. These lie on the neighbouring processes.
* Thus, they have to be copied which is done by this function.
*
* last update 12/02/02, T. Bohlen
*
* ----------------------------------------------------------------------*/
#include "fd.h"
void matcopy_viscac(float ** rho, float ** pi, float ** taup){
extern int MYID, NX, NY, INDEX[5];
extern const int TAG1,TAG2,TAG5,TAG6;
extern FILE *FP;
MPI_Status status;
double time1, time2;
int i, j;
float ** bufferlef_to_rig, ** bufferrig_to_lef;
float ** buffertop_to_bot, ** bufferbot_to_top;
bufferlef_to_rig = matrix(0,NY+1,1,5);
bufferrig_to_lef = matrix(0,NY+1,1,5);
buffertop_to_bot = matrix(0,NX+1,1,5);
bufferbot_to_top = matrix(0,NX+1,1,5);
fprintf(FP,"\n\n **Message from matcopy (written by PE %d):",MYID);
fprintf(FP,"\n Copy material properties at inner boundaries ... \n");
time1=MPI_Wtime();
// if (POS[2]!=0) /* no boundary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
/* storage of top of local volume into buffer */
buffertop_to_bot[i][1] = rho[1][i];
buffertop_to_bot[i][2] = pi[1][i];
buffertop_to_bot[i][5] = taup[1][i];
}
// if (POS[2]!=NPROCY-1) /* no boundary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
/* storage of bottom of local volume into buffer */
bufferbot_to_top[i][1] = rho[NY][i];
bufferbot_to_top[i][2] = pi[NY][i];
bufferbot_to_top[i][5] = taup[NY][i];
}
/*=========sending and receiving of the boundaries==========*/
MPI_Bsend(&buffertop_to_bot[1][1],NX*5,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&buffertop_to_bot[1][1],NX*5,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&status);
MPI_Bsend(&bufferbot_to_top[1][1],NX*5,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferbot_to_top[1][1],NX*5,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&status);
// if (POS[2]!=NPROCY-1) /* no boundary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
rho[NY+1][i] = buffertop_to_bot[i][1];
pi[NY+1][i] = buffertop_to_bot[i][2];
taup[NY+1][i] = buffertop_to_bot[i][5];
}
// if (POS[2]!=0) /* no boundary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
rho[0][i] = bufferbot_to_top[i][1];
pi[0][i] = bufferbot_to_top[i][2];
taup[0][i] = bufferbot_to_top[i][5];
}
// if (POS[1]!=0) /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++)
{
/* storage of left edge of local volume into buffer */
bufferlef_to_rig[j][1] = rho[j][1];
bufferlef_to_rig[j][2] = pi[j][1];
bufferlef_to_rig[j][5] = taup[j][1];
}
// if (POS[1]!=NPROCX-1) /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
/* storage of right edge of local volume into buffer */
bufferrig_to_lef[j][1] = rho[j][NX];
bufferrig_to_lef[j][2] = pi[j][NX];
bufferrig_to_lef[j][5] = taup[j][NX];
}
MPI_Bsend(&bufferlef_to_rig[0][1],(NY+2)*5,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferlef_to_rig[0][1],(NY+2)*5,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&status);
MPI_Bsend(&bufferrig_to_lef[0][1],(NY+2)*5,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferrig_to_lef[0][1],(NY+2)*5,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&status);
// if (POS[1]!=NPROCX-1) /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][NX+1] = bufferlef_to_rig[j][1];
pi[j][NX+1] = bufferlef_to_rig[j][2];
taup[j][NX+1] = bufferlef_to_rig[j][5];
}
// if (POS[1]!=0) /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][0] = bufferrig_to_lef[j][1];
pi[j][0] = bufferrig_to_lef[j][2];
taup[j][0] = bufferrig_to_lef[j][5];
}
if (MYID==0){
time2=MPI_Wtime();
fprintf(FP," finished (real time: %4.2f s).\n",time2-time1);
}
free_matrix(bufferlef_to_rig,0,NY+1,1,5);
free_matrix(bufferrig_to_lef,0,NY+1,1,5);
free_matrix(buffertop_to_bot,0,NX+1,1,5);
free_matrix(bufferbot_to_top,0,NX+1,1,5);
}
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* DENISE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/* $Id: prepare_update_s.c,v 2.2 2007/08/21 13:16:19 tbohlen Exp $*/
#include "fd.h"
void prepare_update_p(float *etajm, float *peta, float **ppi, float **prho, float **ptaup, float **g, float *bjm, float *cjm, float ***e){
extern int NX, NY, L, INVMAT1, MYID;
extern float DT, *FL;
int i, j, l;
extern char MFILE[STRING_SIZE];
extern float F_REF;
float sumpi, ws, *pts;
float pi;
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
}
/*ws=2.0*PI*FL[1];*/
ws=2.0*PI*F_REF;
if(MYID==0)printf("MYID %d: F_REF = %f\n",MYID,F_REF);
sumpi=0.0;
for (l=1;l<=L;l++){
sumpi=sumpi+((ws*ws*pts[l]*pts[l])/(1.0+ws*ws*pts[l]*pts[l]));
}
for (l=1;l<=L;l++){
etajm[l] = peta[l];
}
if (INVMAT1==1){
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
pi=(ppi[j][i]*ppi[j][i]*prho[j][i])/(1.0+sumpi*ptaup[j][i]);
g[j][i] = pi*DT*(1.0+L*ptaup[j][i]);
for (l=1;l<=L;l++){
bjm[l] = 1.0/(1.0+(etajm[l]*0.5));
cjm[l] = 1.0-(etajm[l]*0.5);
e[j][i][l] = pi*etajm[l]*ptaup[j][i];
}
}
}
}
free_vector(pts,1,L);
}
......@@ -387,8 +387,6 @@ void read_par_json(FILE *fp, char *fileinp){
else if(ACOUSTIC==1){
FDORDER=2;
fprintf(fp,"For acoustic modelling only FDORDER=%d possible.\n",FDORDER);
if(L)
err("No viscoacoustic modelling implemented yet");
} else if(ACOUSTIC>1) err("Only ACOUSTIC=0 ((visco-)elastic Modelling and Inversion) or ACOUSTIC=1 (acoustic Modelling and Inversion) possible");
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* DENISE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Read viscoacoustic model properties (vp,density,Qp) from files
*
* Copyright (c) T. Bohlen
* ----------------------------------------------------------------------*/
/* This file contains function readmod, which has the purpose
to read data from model-files for viscoelastic simulation */
#include "fd.h"
void readmod_viscac(float ** rho, float ** pi, float ** taup, float * eta){
extern float DT, *FL, TAU;
extern int L;
extern int NX, NY, NXG, NYG, POS[3], MYID, INVMAT1;
extern char MFILE[STRING_SIZE];
extern FILE *FP;
/* local variables */
float rhov, piv, vp, qp, *pts;
int i, j, ii, jj, l, sw_Qp=1;
FILE *fp_vp, *fp_rho, *fp_qp;
char filename[STRING_SIZE];
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
fprintf(FP,"\n...reading model information from modell-files...\n");
/* read density and seismic velocities */
/* ----------------------------------- */
if(INVMAT1==1){
fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE);
sprintf(filename,"%s.vp",MFILE);
fp_vp=fopen(filename,"r");
if (fp_vp==NULL) err(" Could not open model file for Vp ! ");
fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
sprintf(filename,"%s.rho",MFILE);
fp_rho=fopen(filename,"r");
if (fp_rho==NULL) err(" Could not open model file for densities ! ");
fprintf(FP,"\t Qp:\n\t %s.qp\n\n",MFILE);
sprintf(filename,"%s.qp",MFILE);
fp_qp=fopen(filename,"r");
if (fp_qp==NULL){
if (MYID==0){
printf(" Could not open model file for Qp-values ! \n");
printf(" Uses input file value for tau! \n");
}
sw_Qp=0;
}
}
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
fread(&vp, sizeof(float), 1, fp_vp);
fread(&rhov, sizeof(float), 1, fp_rho);
if (sw_Qp){
fread(&qp, sizeof(float), 1, fp_qp);
}
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
rho[jj][ii]=rhov;
pi[jj][ii]=vp;
if (sw_Qp){
taup[jj][ii]=2.0/qp;}
else taup[jj][ii]=TAU;
}
}
}
fclose(fp_vp);
fclose(fp_rho);
if (sw_Qp) fclose(fp_qp);
/* each PE writes his model to disk */
sprintf(filename,"%s.fdveps.pi",MFILE);
writemod(filename,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(filename,3);
sprintf(filename,"%s.fdveps.rho",MFILE);
writemod(filename,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(filename,3);
free_vector(pts,1,L);
}
......@@ -36,6 +36,7 @@ void surface_acoustic_PML(int ndepth, float ** sp){
j=ndepth; /* The free surface is located exactly in y=1/2*dh !! */
for (i=1;i<=NX;i++){
// sp[j][i] = 0;
for (m=1; m<=fdoh; m++) {
sp[j-m][i] = -sp[j+m][i];
}
......
......@@ -27,13 +27,11 @@
#include "fd.h"
void update_p_PML(int nx1, int nx2, int ny1, int ny2,
float ** vx, float ** vy, float ** sp, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){
void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sp, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){
int i,j, m, fdoh, h, h1;
float g;
float vxx, vyy, vxy, vyx;
......@@ -46,10 +44,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2,
double time1, time2;
dhi = DT/DH;
fdoh = FDORDER/2;
if (infoout && (MYID==0)){
time1=MPI_Wtime();
......@@ -57,64 +53,59 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2,
fprintf(FP," Updating stress components ...");
}