/* * 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, DX, DY, DZ, *FL, TAU; extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID; /*extern char MFILE[STRING_SIZE];*/ char outfile[STRING_SIZE]; /* local variables */ float muv, piv, ws; float Vp, Vs, Rho; float *pts, ts, tp, sumu, sumpi, d,d2,x,y,z, dsteel; int i, j, k, l, ii, jj, kk; /*char modfile[STRING_SIZE];*/ /* background */ const float vp1=5700.0, vs1=3400.0, rho1=2200.0, Qp1=500.0, Qs1=500.0; /* background reflector */ const float vp2=4500.0, vs2=2800.0, rho2=1500.0, Qp2=500.0, Qs2=500.0; const float reflect_dist=80.00; /* tunnel */ const float vpt=0.0, vst=1.0e-4, rhot=1.25, Qpt=1000.0, Qst=1000.0, length_tunnel=60.00, axis_y=(float)NYG*DY/2.0, axis_z=(float)NZG*DZ/2.0, radius=5.0; /*excevation damage zone*/ const float vpedz=1800.0*sqrt(3.0), vsedz=1800.0, rhoedz=1500.0, Qpedz=500.0, Qsedz=500.0, radiusedz=10.0; /* tuebbing */ const float vpconcrete=4000.0, vsconcrete=vpconcrete/sqrt(3), rhoconcrete=2300.0, Qpconcrete=500.0, Qsconcrete=500.0, vpsteel=5860.0, vssteel=3313.0, rhosteel=7700.0, Qpsteel=500.00, Qssteel=500.0, vpbackfill=3000.0, vsbackfill=vpbackfill/sqrt(3), rhobackfill=2000.0, Qpbackfill=500.0, Qsbackfill=500.0, thick_backfill=0.4, thick_tueb=0.6, radius_steel=0.2; /* TBM */ const float vp_tbm_steel=5050, vs_tbm_steel=vp_tbm_steel*sqrt(3), rho_tbm_steel=7700.0, Qp_tbm_steel=500.00, Qs_tbm_steel=500.0, vp_tbm_bentonit=1600.0, vs_tbm_bentonit=1.0e-4, rho_tbm_bentonit=1200.0, Qp_tbm_bentonit=500.00, Qs_tbm_bentonit=500.0, thick_schild=0.3, length_tbm=10.0, thick_cutterhead=0.4, spacing2formation=0.0, thick_bentchamber=2, radius_tbmaxis=0.5; /*-----------------------------------------------------------------------*/ /* 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]; sumu=0.0; sumpi=0.0; 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++){ x=(float)i*DX; y=(float)j*DY; z=(float)k*DZ; /* background */ Vp=vp1; Vs=vs1; Rho=rho1; tp=2.0/Qp1; ts=2.0/Qs1; /* reflector */ if ((dreflect_dist)){ Vp=vp2; Vs=vs2; Rho=rho2; tp=2.0/Qp2; ts=2.0/Qs2;} /*excevation damage zone*/ /*working front*/ d2=sqrt(((y-axis_y)*(y-axis_y))+((z-axis_z)*(z-axis_z))+((x-length_tunnel+radius)*(x-length_tunnel+radius))); if ((d2>radius)&&(d2<(radius+radiusedz))){ Vp=vpedz+((vp1-vpedz)*(d2-radius)/radiusedz); Vs=vsedz+((vs1-vsedz)*(d2-radius)/radiusedz); Rho=rhoedz+((rho1-rhoedz)*(d2-radius)/radiusedz); tp=2.0/Qpedz; ts=2.0/Qsedz;} /*excevation damage zone*/ /*around tube*/ d=sqrt(((y-axis_y)*(y-axis_y))+((z-axis_z)*(z-axis_z))); if ((d>radius)&&(d<(radius+radiusedz))&&(xradius-thick_tueb-thick_backfill)&&(d<(radius-thick_backfill))&&(xradius-thick_backfill)&&(d<(radius))&&(xlength_tunnel-length_tbm)){ Vp=vpt; Vs=vst; Rho=rhot; tp=2.0/Qpt; ts=2.0/Qst;} /*TBM */ /* Schild tube*/ if ((d>radius-thick_schild-spacing2formation)&&(d<(radius-spacing2formation))&&(xlength_tunnel-length_tbm)){ Vp=vp_tbm_steel; Vs=vs_tbm_steel; Rho=rho_tbm_steel; tp=2.0/Qp_tbm_steel; ts=2.0/Qs_tbm_steel;} /*TBM */ /* Schild cover tunnel*/ if ((dlength_tunnel-length_tbm)){ Vp=vp_tbm_steel; Vs=vs_tbm_steel; Rho=rho_tbm_steel; tp=2.0/Qp_tbm_steel; ts=2.0/Qs_tbm_steel;} /*TBM */ /* Schild cover working face*/ if ((dlength_tunnel-thick_bentchamber-thick_schild)){ Vp=vp_tbm_steel; Vs=vs_tbm_steel; Rho=rho_tbm_steel; tp=2.0/Qp_tbm_steel; ts=2.0/Qs_tbm_steel;} /*TBM */ /* cutter head */ if ((dlength_tunnel-thick_cutterhead-spacing2formation)){ Vp=vp_tbm_steel; Vs=vs_tbm_steel; Rho=rho_tbm_steel; tp=2.0/Qp_tbm_steel; ts=2.0/Qs_tbm_steel;} /*TBM */ /* axis connecting tbm schild body and cutter head */ if ((dlength_tunnel-thick_cutterhead-spacing2formation-thick_bentchamber)){ Vp=vp_tbm_steel; Vs=vs_tbm_steel; Rho=rho_tbm_steel; tp=2.0/Qp_tbm_steel; ts=2.0/Qs_tbm_steel;} 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; 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; } } } } sprintf(outfile,"model/tuebbingtbm.vs"); writemod(outfile,u,3); MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(outfile,3); sprintf(outfile,"model/tuebbingtbm.vp"); writemod(outfile,pi,3); MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(outfile,3); /*sprintf(outfile,"model/test.u"); writemod(outfile,u,3); MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(outfile,3); */ sprintf(outfile,"model/tuebbingtbm.rho"); writemod(outfile,rho,3); MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(outfile,3); /* sprintf(outfile,"model/test.pi"); writemod(outfile,pi,3); MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(outfile,3); */ free_vector(pts,1,L); }