new toy example script

script runs forward and inversion

par/in_and_out/ifos3d_toy.json

@@ -32,7 +32,6 @@
 "Model" : "comment",
 			"READMOD" : "0",
 			"MFILE" : "model/toy",
-			
 
 "Q-approximation" : "comment",
 			"L" : "0",
@@ -60,19 +59,15 @@
 			"REC_ARRAY_DIST" : "30.0",
 			"DRX" : "10",
 			"DRY" : "10",
-			
 
 "Seismograms" : "comment",
 			"NDT" : "1",
 			"SEIS_FORMAT" : "1",
 			"SEIS_FILE" : "./su/cal_toy",
 
-			"VERBOSE" : "0",
-	
-	
 "Method" : "comment",
-			"METHOD" : "0",
-			
+			"METHOD" : "1",
+
 "INVERSION PARAMETERS" : "comment",
 
 "In- and Output Files" : "comment",

par/run_toy_example.sh

+# run forward simulation to obtain observed data
+make clean
+make install  MODEL=hh_toy_true.c
+mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/ifos3d_toy_FW.json | tee ./in_and_out/ifos3D.out
+
+# invert observed data using homogeneous starting model
+make clean
+make install MODEL=hh_toy_start.c
+ mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/ifos3d_toy.json | tee ./in_and_out/ifos3D_INV.out

src/Makefile

@@ -8,7 +8,7 @@
 
 # source code for model generation (elastic)
 
-MODEL_SCR = hh_toy.c
+MODEL_SCR = hh_toy_true.c
 
 # Compiler (LAM: CC=hcc, CRAY T3E: CC=cc, SHARCNet: mpicc)
 

src/hh_toy_start.c

+/*
+ *   homogeneous half space 
+ *   last update 02.11.02, T. Bohlen
+ */
+
+#include "fd.h"
+
+void model(float  ***  rho, float ***  pi, float ***  u, 
+float ***  taus, float ***  taup, float *  eta){
+
+	/*-----------------------------------------------------------------------*/
+	/* extern variables */
+
+	extern float DT, *FL, TAU;
+	extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L;
+	extern FILE *FP;
+	/* local variables */
+	float muv, piv, ws;
+	float Vp, Vs, Rho;
+	float *pts=NULL, ts=0.0, tp=0.0, sumu, sumpi;
+	int i, j, k, l, ii, jj, kk;
+	int kasten=0;
+
+	/* parameters for layer 1 */
+	const float vp1=6200.0, vs1=3600.0, rho1=2800.0;
+	/*const float vp2=6200.0, vs2=3600.0, rho2=2800.0;*/
+      /*const float vp2=7000.0, vs2=3900.0, rho2=2800.0;*/
+
+	/*-----------------------------------------------------------------------*/
+	sumu=0.0; 
+	sumpi=0.0;
+
+	fprintf(FP," start creation of toy-model");
+	
+	 if(L){
+		/* 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];
+		}
+		ts=TAU;  
+		tp=TAU;
+		ws=2.0*PI*FL[1];
+		for (l=1;l<=L;l++){
+			sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
+			sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
+		}
+	 }
+
+	/* loop over global grid */
+	for (k=1;k<=NZG;k++){
+		for (i=1;i<=NXG;i++){
+			for (j=1;j<=NYG;j++){
+
+				
+				Vp=vp1; Vs=vs1; Rho=rho1; 
+				
+				/*for usability reasons, "z" - as commonly used - denotes the depth (vertical direction),
+				  however, internally "y" is used for the vertical coordinate,
+				  we simply switch the "y" and "z" coordinate as read in the input file,
+				  therefore use the variable "j" to, e.g., calculate/address the vertical coordinate:
+				  x=(float)i*DX;
+				  y=(float)k*DZ;
+				  z=(float)j*DY;*/  				
+
+				muv=Vs*Vs*Rho/(1.0+sumu);
+				piv=Vp*Vp*Rho/(1.0+sumpi);
+
+				/* 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)) && 
+				    (POS[3]==((k-1)/NZ))){
+					ii=i-POS[1]*NX;
+					jj=j-POS[2]*NY;
+					kk=k-POS[3]*NZ;
+					if(L){
+						taus[jj][ii][kk]=ts;
+						taup[jj][ii][kk]=tp;
+					}
+					u[jj][ii][kk]=muv;
+					rho[jj][ii][kk]=Rho;
+					pi[jj][ii][kk]=piv;
+				}
+			}
+		}
+	}
+if(kasten){
+	for (k=50;k<=90;k++){
+		for (i=50;i<=80;i++){
+			for (j=50;j<=65;j++){	/*vertical*/
+
+				Vp=0.0; Vs=0.0; Rho=0.0;
+				Vp=vp1+500.0; Vs=vs1+300; Rho=rho1; 
+
+				muv=Vs*Vs*Rho/(1.0+sumu);
+				piv=Vp*Vp*Rho/(1.0+sumpi);
+
+				/* 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)) && 
+				    (POS[3]==((k-1)/NZ))){
+					ii=i-POS[1]*NX;
+					jj=j-POS[2]*NY;
+					kk=k-POS[3]*NZ;
+					if(L){
+						taus[jj][ii][kk]=ts;
+						taup[jj][ii][kk]=tp;
+					}
+					u[jj][ii][kk]=muv;
+					rho[jj][ii][kk]=Rho;
+					pi[jj][ii][kk]=piv;
+				}
+			}
+			
+		}
+	}
+
+	for (k=50;k<=90;k++){
+		for (i=81;i<=110;i++){
+			for (j=50;j<=95;j++){
+
+				Vp=0.0; Vs=0.0; Rho=0.0;
+				Vp=vp1-250; Vs=vs1-100; Rho=rho1; 
+
+				muv=Vs*Vs*Rho/(1.0+sumu);
+				piv=Vp*Vp*Rho/(1.0+sumpi);
+
+				/* 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)) && 
+				    (POS[3]==((k-1)/NZ))){
+					ii=i-POS[1]*NX;
+					jj=j-POS[2]*NY;
+					kk=k-POS[3]*NZ;
+					if(L){
+						taus[jj][ii][kk]=ts;
+						taup[jj][ii][kk]=tp;
+					}
+					u[jj][ii][kk]=muv;
+					rho[jj][ii][kk]=Rho;
+					pi[jj][ii][kk]=piv;
+				}
+			}
+			
+		}
+	}	for (k=50;k<=90;k++){
+		for (i=50;i<=80;i++){
+			for (j=66;j<=95;j++){
+
+				Vp=0.0; Vs=0.0; Rho=0.0;
+				Vp=vp1+300; Vs=vs1+200; Rho=rho1; 
+
+				muv=Vs*Vs*Rho/(1.0+sumu);
+				piv=Vp*Vp*Rho/(1.0+sumpi);
+
+				/* 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)) && 
+				    (POS[3]==((k-1)/NZ))){
+					ii=i-POS[1]*NX;
+					jj=j-POS[2]*NY;
+					kk=k-POS[3]*NZ;
+					if(L){
+						taus[jj][ii][kk]=ts;
+						taup[jj][ii][kk]=tp;
+					}
+					u[jj][ii][kk]=muv;
+					rho[jj][ii][kk]=Rho;
+					pi[jj][ii][kk]=piv;
+				}
+			}
+			
+		}
+	}	for (k=91;k<=110;k++){
+		for (i=50;i<=110;i++){
+			for (j=50;j<=95;j++){
+
+				Vp=0.0; Vs=0.0; Rho=0.0;
+				Vp=vp1-400; Vs=vs1-150; Rho=rho1;  				
+
+				muv=Vs*Vs*Rho/(1.0+sumu);
+				piv=Vp*Vp*Rho/(1.0+sumpi);
+
+				/* 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)) && 
+				    (POS[3]==((k-1)/NZ))){
+					ii=i-POS[1]*NX;
+					jj=j-POS[2]*NY;
+					kk=k-POS[3]*NZ;
+					if(L){
+						taus[jj][ii][kk]=ts;
+						taup[jj][ii][kk]=tp;
+					}
+					u[jj][ii][kk]=muv;
+					rho[jj][ii][kk]=Rho;
+					pi[jj][ii][kk]=piv;
+				}
+			}
+			
+		}
+	}
+	
+}
+	
+	
+	MPI_Barrier(MPI_COMM_WORLD);
+         if(L)	free_vector(pts,1,L);
+}
+
+
+