Commit faee198a authored by niklas.thiel's avatar niklas.thiel

initial double-couple implementation

parent 15c13ee0
......@@ -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",
......
#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;
......
......@@ -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 ! ");
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment