...
 
Commits (6)
......@@ -99,7 +99,7 @@ The program checks these criteria for the entire model, outputs a warning messag
"SOURCE_TYPE" : "3",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;
rotated_force=4" : "comment",
rotated_force=4;double-couple=5" : "comment",
"SRCREC" : "1",
"SRCREC values : read source positions from SOURCE_FILE=1,
......@@ -187,6 +187,8 @@ The wavelets in each su file must have the same number of samples as specified b
The following source types are availabe: explosive sources that excite compressional waves only (SOURCE\_TYPE=1), and point forces in the x- and y-direction (SOURCE\_TYPE=2,3).
The force sources excite both P- and S-waves. The explosive source is located at the same position as the diagonal elements of the stress tensor, i.e. at (i,j) (Figure \ref{fig_cell}).
The forces are located at the same position as the corresponding components of particle velocity (Figure \ref{fig_cell}). If (x,y) denotes the position at which the source location is defined in source.dat, then the actual force in x-direction is located at (x+DX/2,y) and the actual force in y-direction is located at (x,y+DY/2). With SOURCE\_TYPE=4 a custom directive force can be defined by a force angle between y and x. The angle of the force must be specified in the SOURCE\_FILE after AMP. This force is not aligned along the main directions.
By using SOURCE\_TYPE=5, a 2D double-couple source is activated. The $M_0$ of the moment tensor is given in AMP in the SOURCE\_FILE, the angle in degree also in the SOURCE\_FILE at the SOURCE\_AZIMUTH position.
\newline
The locations of multiple sources must be defined in an external ASCII file (SOURCE\_FILE) that has the following format:
......@@ -199,7 +201,7 @@ In the following lines, you can define certain parameters for each source point:
The first line must be the overall number of sources (NSRC). XSRC is the x-coordinate of a source point (in meter), YSRC is the y-coordinate of a source point (in meter). ZSRC is the z-coordinate should always be set to 0.0, because IFOS2D is a 2D code. TD is the excitation time (time-delay) for the source point (in seconds), FC is the center frequency of the source signal (in Hz), and AMP is the maximum amplitude of the source signal.
\newline
\textbf{Optional parameter:} The SOURCE\_AZIMUTH if SOURCE\_TYPE is 4. The SOURCE\_AZIMUTH is the angle between the y- and x-direction in degree and with SOURCE\_TYPE if SOURCE\_TYPE is set here, the value of SOURCE\_TYPE in the input file is ignored.
\textbf{Optional parameter:} The SOURCE\_AZIMUTH if SOURCE\_TYPE is 4 or 5. The SOURCE\_AZIMUTH is the angle between the y- and x-direction in degree and with SOURCE\_TYPE if SOURCE\_TYPE is set here, the value of SOURCE\_TYPE in the input file is ignored.
The SOURCE\_FILE = ./sources/source.dat that defines an explosive source at $x_s=2592.0\;$ m and $y_s=2106.0\;$ m with
a center frequency of 5 Hz (no time delay) is
......@@ -489,8 +491,11 @@ Default values are:
With the use of a workflow file, you can define different FWI stages. For instance, one FWI stage could refer to one corner frequency of a low-pass filter or to a specific time window. Every line in the workflow file corresponds to one FWI stage. To use a workflow file the switch USE\_WORKFLOW have to be set to 1. The algorithm will automatically change to the next line of the workflow file if the abort criterium of the current line is reached or if no step length could be found which reduces the misfit. The structure of the variables inside the workflow file is as follow: \\
\noindent
\begin{tabular}{llllllllllIlll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE & GAMMA \\
\begin{tabular}{lllllllll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE\\
\end{tabular}
\begin{tabular}{llllll}
\# & 0 & 0 & EPRECOND & EPSILON\_WE & GAMMA \\
\end{tabular}
\ \\
......
......@@ -8,13 +8,13 @@
"MAXRELERROR" : "0",
"2-D Grid" : "comment",
"NX" : "500",
"NY" : "100",
"DH" : "0.2",
"NX" : "4800",
"NY" : "2040",
"DH" : "50",
"Time Stepping" : "comment",
"TIME" : "0.5",
"DT" : "5.0e-05",
"TIME" : "20",
"DT" : "4.0e-03",
"Source" : "comment",
"SOURCE_SHAPE" : "4",
......@@ -22,31 +22,31 @@
"SOURCE_SHAPE values: Spike=6;from_SIGNAL_FILE_in_su_format=7;integral_of_SIN**3=8" : "comment",
"SIGNAL_FILE" : "./ormsby.dat",
"SOURCE_TYPE" : "3",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4" : "comment",
"SOURCE_TYPE" : "5",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4;double-couple=5" : "comment",
"SOURCE_FILE" : "./source/sources.dat",
"SOURCE_FILE" : "./source/sources_dc1.dat",
"RUN_MULTIPLE_SHOTS" : "1",
"Model" : "comment",
"READMOD" : "0",
"MFILE" : "model/model_Test",
"READMOD" : "1",
"MFILE" : "model_true/ecusub_2D_0N_4800_2040_50",
"Free Surface" : "comment",
"FREE_SURF" : "1",
"PML Boundary" : "comment",
"FW" : "20",
"VPPML" : "600.0",
"FPML" : "31.25",
"FW" : "15",
"VPPML" : "5000.0",
"FPML" : "2",
"BOUNDARY" : "0",
"npower" : "4.0",
"k_max_PML" : "8.0",
"Receiver" : "comment",
"SEISMO" : "1",
"SEISMO" : "5",
"READREC" : "1",
"REC_FILE" : "./receiver/receiver.dat",
"REC_FILE" : "./receiver/rec_eqsource.dat",
"REFRECX, REFRECY" : "0.0 , 0.0",
"XREC1, YREC1" : "6.0 , 0.2",
"XREC2, YREC2" : "93.0 , 0.2",
......@@ -56,6 +56,17 @@
"SEIS_FORMAT" : "1",
"SEIS_FILE" : "su/IFOS",
"Snapshots" : "comment",
"SNAP" : "4",
"TSNAP1" : "4e-3",
"TSNAP2" : "12.0",
"TSNAPINC" : "0.12",
"IDX" : "1",
"IDY" : "1",
"SNAP_FORMAT" : "3",
"SNAPSHOT_START , SNAPSHOT_END , SNAPSHOT_INCR" : "1, 1 , 1",
"SNAP_FILE" : "./snap/waveform_forward",
"General inversion parameters" : "comment",
"PARAMETERIZATION" : "1",
"FORWARD_ONLY" : "10",
......
......@@ -23,7 +23,7 @@
"SIGNAL_FILE" : "./ormsby.dat",
"SOURCE_TYPE" : "3",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4" : "comment",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4;double-couple=5" : "comment",
"SOURCE_FILE" : "./source/sources.dat",
"RUN_MULTIPLE_SHOTS" : "1",
......
......@@ -23,7 +23,7 @@
"SIGNAL_FILE" : "./ormsby.dat",
"SOURCE_TYPE" : "3",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4" : "comment",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4;double-couple=5" : "comment",
"SOURCE_FILE" : "./source/sources.dat",
"RUN_MULTIPLE_SHOTS" : "1",
......
......@@ -23,7 +23,7 @@
"SIGNAL_FILE" : "./STF/stf.su",
"SOURCE_TYPE" : "3",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4" : "comment",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4;double-couple=5" : "comment",
"SOURCE_FILE" : "./source/sources.dat",
"RUN_MULTIPLE_SHOTS" : "1",
......
#FLNODE
#converted by MOCOX V1.2 from final.mod.gremlin
# 0.0000 1 0 6371.0000
# depth rho vp vs
# in m in g/cm^3 in km/s in km/s
#7
0.0000 1.7030 0.3569 0.0917
0.9299 1.7030 0.3569 0.1661
0.9299 1.6960 0.3569 0.1661
6.2636 1.6960 0.3569 0.2479
6.2636 2.0010 1.7450 0.2735
6.3800 2.0010 1.7450 0.2751
6.3800 2.0010 1.7450 0.2751
#FLNODE
#converted by MOCOX V1.2 from final.mod.gremlin
# 0.0000 1 0 6371.0000
# depth rho vp vs
# in m in g/cm^3 in km/s in km/s
#3
0.0000 1.7030 0.3569 0.0917
9.0000 2.0010 1.7450 0.2751
9.0000 2.0010 1.7450 0.2751
#FLNODE
#converted by MOCOX V1.2 from final.mod.gremlin
# 0.0000 1 0 6371.0000
# depth rho vp
# in m in g/cm^3 in km/s
#6
0.0000 1.0200 1.4800
500.0000 1.0200 1.4800
500.0000 2.0010 1.5000
1000.0000 2.0010 2.0000
1000.0000 2.0010 2.1000
2000.0000 2.0010 2.6000
#FLNODE
#converted by MOCOX V1.2 from final.mod.gremlin
# 0.0000 1 0 6371.0000
# depth rho vp taper
# in m in g/cm^3 in km/s
#5
0.0000 1.0200 1.4800 0
500.0000 1.0200 1.4800 0
500.0000 2.0010 1.5000 0
550.0000 2.0010 1.5000 1
2000.0000 2.0010 2.5000 1
rm -f su/IFOS* su/measured_data/IFOS*
#rm -f su/IFOS* su/measured_data/IFOS*
lamboot
#lamboot
#lamboot -v lamhosts
mpirun -np 8 nice -19 ../bin/IFOS2D in_and_out/IFOS_FW.json | tee in_and_out/IFOS_FW.out
cd model
cp model_Test_rho_it_0.bin model_Test_vp_it_0.bin model_Test_vs_it_0.bin ../model_true/.
cd ..
mpirun -np 8 nice -19 ../bin/IFOS2D in_and_out/IFOS2D_FW.json | tee in_and_out/IFOS2D_FW.out
#cd model
#cp model_Test_rho_it_0.bin model_Test_vp_it_0.bin model_Test_vs_it_0.bin ../model_true/.
#cd ..
cd su
#cd su
for ((i=1; i < ($1+1); i++)) ; do
#for ((i=1; i < ($1+1); i++)) ; do
mv IFOS_y.su.shot$i.it1 measured_data/IFOS_y.su.shot$i
mv IFOS_x.su.shot$i.it1 measured_data/IFOS_x.su.shot$i
# mv IFOS_y.su.shot$i.it1 measured_data/IFOS_y.su.shot$i
# mv IFOS_x.su.shot$i.it1 measured_data/IFOS_x.su.shot$i
done
#done
cd ..
#cd ..
......@@ -1269,7 +1269,10 @@ int main(int argc, char **argv){
/* explosive source */
if ((SOURCE_TYPE==1))
psource(nt,psxx,psyy,psp,srcpos_loc,signals,nsrc_loc,0);
/* earthquake source */
if ((SOURCE_TYPE==5))
eqsource(nt,psxx,psyy,psxy,psp,srcpos_loc,signals,nsrc_loc,0);
/* Applying free surface condition */
if ((FREE_SURF) && (POS[2]==0)){
......@@ -1713,6 +1716,10 @@ int main(int argc, char **argv){
/* explosive source */
if ((SOURCE_TYPE==1)&&(WAVETYPE==1||WAVETYPE==3))
psource(nt,psxx,psyy,psp,srcpos_loc,signals,nsrc_loc,0);
/* earthquake source */
if ((SOURCE_TYPE==5))
eqsource(nt,psxx,psyy,psxy,psp,srcpos_loc,signals,nsrc_loc,0);
/* Applying free surface condition */
if ((FREE_SURF) && (POS[2]==0)){
......@@ -3462,6 +3469,10 @@ int main(int argc, char **argv){
/* explosive source */
if ((SOURCE_TYPE==1))
psource(nt,psxx,psyy,psp,srcpos_loc,signals,nsrc_loc,0);
/* earthquake source */
if ((SOURCE_TYPE==5))
eqsource(nt,psxx,psyy,psxy,psp,srcpos_loc,signals,nsrc_loc,0);
/* Applying free surface condition */
if ((FREE_SURF) && (POS[2]==0)){
......
......@@ -84,6 +84,7 @@ IFOS2D= \
psource.c \
holbergcoeff.c\
comm_ini.c\
eqsource.c\
exchange_v.c \
exchange_s.c \
exchange_p.c \
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2016 For the list of authors, see file AUTHORS.
*
* This file is part of IFOS.
*
* IFOS 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.
*
* IFOS 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 IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Source excitation using Moment tensor for simulation of earthquakes (double-couple)
* ----------------------------------------------------------------------*/
#include "fd.h"
void eqsource(int nt, float ** sxx, float ** syy, float ** sxy, float ** sp,
float ** srcpos_loc, float ** signals, int nsrc, int sw){
extern int MYID, NT, ACOUSTIC;
extern float DT,DH;
int i, j, l;
float amp,dip,amon,dip1,scale_amp;
float m11,m22,m12;
/* adding source wavelet to stress components
(moment tensor source) at source points */
for (l=1;l<=nsrc;l++) {
i=(int)srcpos_loc[1][l];
j=(int)srcpos_loc[2][l];
amon=srcpos_loc[6][l];
dip1=srcpos_loc[7][l];
dip=dip1*PI/180.0;
scale_amp=DT/(DH*DH);
amp = amon*scale_amp*signals[l][nt];
m11 = - sin(2*dip);
m12 = - cos(2*dip);
m22 = sin(2*dip);
sxx[j][i] += amp * m11;
sxy[j][i] += 0.3 * amp * m12;
sxy[j-1][i] += 0.3 * amp * m12;
sxy[j][i-1] += 0.3 * amp * m12;
syy[j][i] += amp * m22;
}
}
......@@ -103,6 +103,9 @@ void dealloc_sections(int ntr,int ns,int **recpos_loc,float **sectionvx,float **
float **sectionpdata,float **sectionpdiff,float **sectionpdiffold,float **sectionvxdata,float **sectionvxdiff,float **sectionvxdiffold,float **sectionvydata,
float **sectionvydiff,float **sectionvydiffold,float **sectionvzdata,float **sectionvzdiff,float **sectionvzdiffold);
void eqsource(int nt, float ** sxx, float ** syy, float ** sxy, float ** sp,
float ** srcpos_loc, float ** signals, int nsrc, int sw);
float exchange_L2(float L2, int sw, int bcast_l2);
void exchange_rsg(float ** vx, float ** vy, float ** vz,
......
......@@ -189,7 +189,7 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("SOURCE_TYPE",number_readobjects,&SOURCE_TYPE,varname_list, value_list))
declare_error("Variable SOURCE_TYPE could not be retrieved from the json input file!");
/* Definition of inversion for source time function */
if (get_int_from_objectlist("INV_STF",number_readobjects,&INV_STF,varname_list, value_list)){
INV_STF=0;
......
......@@ -143,17 +143,8 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
}
if(FORWARD_ONLY==1){
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
}
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.shot%d",SEIS_FILE,ishot);
}
......@@ -165,7 +156,7 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
switch (SEISMO){
case 1 : /* particle velocities only */
if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(ADJOINT_TYPE==1 || ADJOINT_TYPE==0) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
......@@ -201,7 +192,7 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 4 : /* everything */
if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(ADJOINT_TYPE==1 || ADJOINT_TYPE==0) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
......@@ -232,7 +223,7 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 5 : /* everything except curl and div */
if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(ADJOINT_TYPE==1 || ADJOINT_TYPE==0) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
......
......@@ -81,10 +81,23 @@ int main(int argc, char **argv){
if(WAVETYPE==2||WAVETYPE==3) {
merge(nsnap,7);
}
/*pressure */
merge(nsnap,6);
case 3 :
merge(nsnap,4);
merge(nsnap,5);
break;
case 5 :
if(WAVETYPE==1||WAVETYPE==3) {
merge(nsnap,1);
merge(nsnap,2);
}
if(WAVETYPE==2||WAVETYPE==3) {
merge(nsnap,7);
}
/*pressure */
merge(nsnap,6);
break;
default :
warning(" snapmerge: cannot identify content of snapshot !");
break;
......
This diff is collapsed.
......@@ -120,8 +120,8 @@ void write_par(FILE *fp){
fprintf(fp," ------------------------- SOURCE -----------------------------\n");
if (SRCREC){
fprintf(fp," reading source positions, time delay, centre frequency \n");
fprintf(fp," and initial amplitude from ASCII-file \n");
fprintf(fp," reading source positions, time delay, centre frequency, \n");
fprintf(fp," initial amplitude and source azimuth from ASCII-file \n");
fprintf(fp,"\t%s\n\n",SOURCE_FILE);
} else {
fprintf(fp," plane wave excitation: depth= %5.2f meter \n",PLANE_WAVE_DEPTH);
......@@ -176,6 +176,9 @@ void write_par(FILE *fp){
case 4 :
fprintf(fp," rotated point source with directive force in x- and y-direction\n");
break;
case 5 :
fprintf(fp," double-couple source\n");
break;
default :
declare_error(" Sorry, wrong source type specification ! ");
}
......