diff --git a/src/EW.h b/src/EW.h index f59ac963..6d2037cf 100644 --- a/src/EW.h +++ b/src/EW.h @@ -119,7 +119,7 @@ void solve( vector & a_GlobalSources, vector & a_GlobalTim vector& U, vector& Um, vector& Upred_saved_sides, vector& Ucorr_saved_sides, bool save_sides, int event, int save_steps, - int varcase, vector& pseudoHessian ); + int varcase, vector& pseudoHessian, int step_to_record); void solveTT( Source* a_GlobalSource, vector & a_GlobalTimeSeries, double* xs, int nmpars, MaterialParameterization* mp, int wave_mode, int event, int myrank); @@ -134,7 +134,7 @@ void solve_backward_allpars( vector & a_GlobalSources, vector& vector& a_Lambda, vector & a_GlobalTimeSeries, vector& a_U, vector& a_Um, vector& Upred_saved_sides, vector& Ucorr_saved_sides, float_sw4 gradients[11], - vector& gRho, vector& gMu, vector& gLambda, int event ); + vector& gRho, vector& gMu, vector& gLambda, int step_to_record, int event ); //int nmpar, float_sw4* gradientm ); void solve_dudp( vector& a_Sources, vector& a_Rho, diff --git a/src/Mopt.C b/src/Mopt.C index 289ee620..19ff8420 100644 --- a/src/Mopt.C +++ b/src/Mopt.C @@ -55,6 +55,7 @@ Mopt::Mopt( EW* a_ew ) m_vs_min = -100.; m_vs_max = -100.; m_freq_peakpower=0.0; + m_skip_precursor=false; m_wave_mode=2; // default to both P and S waves otherwise 0 for P and 1 for S only m_win_mode =0; // default, use windows stored on hdf5-file. m_twin_shift=-0.5; @@ -534,6 +535,13 @@ void Mopt::processMrun( char* buffer ) } m_ew->set_filterpar(filterpar); } + else if( startswith("skip_precursor=",token) ) + { + token += 15; + int n = strlen(token); + m_skip_precursor = strncmp("yes",token,n)== 0 || strncmp("true",token,n)==0 || strncmp("1",token,n)== 0; + if(m_myrank == 0) cout << "skip_precursor=" << m_skip_precursor << endl; + } else if( startswith("wave_mode=",token) ) { token += 10; diff --git a/src/Mopt.h b/src/Mopt.h index 4d340646..b17a6d91 100644 --- a/src/Mopt.h +++ b/src/Mopt.h @@ -54,6 +54,7 @@ class Mopt float_sw4 get_vs_min() const { return m_vs_min; } float_sw4 get_vs_max() const { return m_vs_max; } float_sw4 get_freq_peakpower() const { return m_freq_peakpower; } + bool get_skip_precursor() const { return m_skip_precursor; } int get_wave_mode() const { return m_wave_mode; } float get_twin_shift() const { return m_twin_shift; } float get_twin_scale() const { return m_twin_scale; } @@ -80,6 +81,7 @@ class Mopt // FWI workflow options float_sw4 m_vp_min, m_vp_max, m_vs_min, m_vs_max; // global velocity constraints float_sw4 m_freq_peakpower; // peak-power freq for setting traveltime windows or smoothing gradients + bool m_skip_precursor; // skip BC recording for source precursor to save compute time and memory int m_wave_mode; // 0: P 1: S 2: both float m_twin_shift, m_twin_scale; int m_win_mode; diff --git a/src/Source.h b/src/Source.h index 766ac34d..c24184a9 100644 --- a/src/Source.h +++ b/src/Source.h @@ -63,7 +63,7 @@ class Source const char *name, bool topodepth, int ncyc=1, - float_sw4* pars=NULL, int npars=0, int* ipars=NULL, int nipars=0, bool correctForMu=false ); + float_sw4* pars=NULL, int npars=0, int* ipars=NULL, int nipars=0, bool correctForMu=false); Source(EW * a_ew, float_sw4 frequency, float_sw4 t0, float_sw4 x0, float_sw4 y0, float_sw4 z0, @@ -74,7 +74,7 @@ class Source const char *name, bool topodepth, int ncyc=1, - float_sw4* pars=NULL, int npars=0, int* ipars=NULL, int nipars=0, bool correctForMu=false ); + float_sw4* pars=NULL, int npars=0, int* ipars=NULL, int nipars=0, bool correctForMu=false); ~Source(); @@ -134,7 +134,7 @@ class Source bool get_CorrectForMu(){return mShearModulusFactor;}; void set_CorrectForMu(bool smf){mShearModulusFactor=smf;}; float_sw4 getTimeOffset() const { return mT0; }; - + private: Source(); diff --git a/src/TimeSeries.h b/src/TimeSeries.h index 19e28b45..3c21692d 100644 --- a/src/TimeSeries.h +++ b/src/TimeSeries.h @@ -265,7 +265,6 @@ float_sw4 m_scalefactor; // Window for optimization, m_winL, m_winR given relative simulation time zero. float_sw4 m_winL, m_winR, m_winL2, m_winR2; bool m_use_win, m_use_x, m_use_y, m_use_z; - // quiet mode? bool mQuietMode; diff --git a/src/lbfgs.C b/src/lbfgs.C index 56d7d04a..444e4e10 100644 --- a/src/lbfgs.C +++ b/src/lbfgs.C @@ -621,7 +621,8 @@ void lbfgs( EW& simulation, int nspar, int nmpars, double* xs, { #ifdef USE_HDF5 // Tang: need to create a HDF5 file before writing - if (GlobalTimeSeries[e].size() > 0 && GlobalTimeSeries[e][0]->getUseHDF5()) { + if (GlobalTimeSeries[e].size() > 0 && GlobalTimeSeries[e][0]->getUseHDF5()) + { for (int tsi = 0; tsi < GlobalTimeSeries[e].size(); tsi++) GlobalTimeSeries[e][tsi]->resetHDF5file(); if( GlobalTimeSeries[e][0]->myPoint() ) @@ -633,6 +634,7 @@ void lbfgs( EW& simulation, int nspar, int nmpars, double* xs, #endif for( int m = 0; m < GlobalTimeSeries[e].size(); m++ ) GlobalTimeSeries[e][m]->writeFile( "_ini" ); + } } if( myRank == 0 ) @@ -856,6 +858,12 @@ void lbfgs( EW& simulation, int nspar, int nmpars, double* xs, dam[i] = dm[i]; int retcode; double fp; + + // save time windows per iteration + for( int e=0 ; e < GlobalTimeSeries.size(); e++ ) + for( int m=0 ; m < GlobalTimeSeries[e].size() ; m++ ) + GlobalTimeSeries[e][m]->writeWindows(); + if( myRank == 0 ) cout << "Line search.. " << endl; @@ -1052,6 +1060,8 @@ void lbfgs( EW& simulation, int nspar, int nmpars, double* xs, dfs[i] = dfps[i]; for( int i=0 ; i < nmpard ; i++ ) dfm[i] = dfpm[i]; + + } // end while it U(ng), Um(ng), ph(ng); simulation.solve( GlobalSources[0], GlobalTimeSeries[0], simulation.mMu, simulation.mLambda, simulation.mRho, U, Um, upred_saved, - ucorr_saved, false, 0, 0, 0, ph ); + ucorr_saved, false, 0, 0, 0, ph, 1); // save all time series diff --git a/src/moptmain.C b/src/moptmain.C index 677671bb..7985da00 100644 --- a/src/moptmain.C +++ b/src/moptmain.C @@ -225,6 +225,33 @@ void set_timewindows_from_eikonal_time(vector >& GlobalTimeS } // end of events } +int firstTimeStep_to_recordBC(const vector & a_Sources, const vector & a_TimeSeries, const Mopt *mopt, const float_sw4 dt, +const int e, const int myrank) +{ + if(!mopt->get_skip_precursor()) return 1; + + //float_sw4 fpeak = mopt->get_freq_peakpower() > 0? mopt->get_freq_peakpower() : a_Sources[0]->getFrequency(); + //int step_to_record = a_Sources[0]->getTimeOffset()/dt- floor(1./fpeak/dt*1.5); + + FILE *ft; + float tp, ts, t_trunc; + char file[STRINGSIZE]; + sprintf(file, "%s/time_to_record_%d.txt", a_TimeSeries[0]->getPath().c_str(),e); + ft = fopen(file, "r"); + if(ft!=NULL) { + fscanf(ft, "%g\t%g\n", &tp, &ts); + if(myrank==0) printf("tp=%g\tts=%g\n", tp, ts); + } + + t_trunc = mopt->get_wave_mode()==1? ts : tp; + + //time to record BC is based on the minimum time to reach BC's if provided in traveltime preprocessing or a default fraction of source t0 + int step_to_record = ft==NULL? int(a_Sources[0]->getTimeOffset()/dt*0.38) : int(t_trunc*0.9/dt); + if(myrank==0) printf("t_trunc=%g\tstep-to-record=%d out of two truncation criteria %d and %d\n", t_trunc, step_to_record, + int(a_Sources[0]->getTimeOffset()/dt*0.38), int(t_trunc*0.9/dt)); + return step_to_record>1 ? step_to_record : 1; +} + //----------------------------------------------------------------------- void compute_f( EW& simulation, int nspar, int nmpars, double* xs, int nmpard, double* xm, @@ -334,7 +361,10 @@ void compute_f( EW& simulation, int nspar, int nmpars, double* xs, { // simulation.solve( src, GlobalTimeSeries[e], mu, lambda, rho, U, Um, upred_saved, ucorr_saved, false, e ); sw4_profile->time_stamp("forward solve" ); - simulation.solve( GlobalSources[e], GlobalTimeSeries[e], mu, lambda, rho, U, Um, upred_saved, ucorr_saved, false, e, mopt->m_nsteps_in_memory, 0, ph ); + int myrank; + MPI_Comm_rank(MPI_COMM_WORLD,&myrank); + int step_to_record = firstTimeStep_to_recordBC(GlobalSources[e], GlobalTimeSeries[e], mopt, simulation.getTimeStep(), e, myrank); + simulation.solve( GlobalSources[e], GlobalTimeSeries[e], mu, lambda, rho, U, Um, upred_saved, ucorr_saved, false, e, mopt->m_nsteps_in_memory, 0, ph, step_to_record); sw4_profile->time_stamp("done forward solve" ); // cout.precision(16); // Compute misfit @@ -564,7 +594,8 @@ void compute_f_and_df( EW& simulation, int nspar, int nmpars, double* xs, for( int e=0 ; e < simulation.getNumberOfLocalEvents() ; e++ ) { sw4_profile->time_stamp("forward solve" ); - simulation.solve( GlobalSources[e], GlobalTimeSeries[e], mu, lambda, rho, U, Um, upred_saved, ucorr_saved, true, e, mopt->m_nsteps_in_memory, phcase, pseudo_hessian ); + int step_to_record = firstTimeStep_to_recordBC(GlobalSources[e], GlobalTimeSeries[e], mopt, simulation.getTimeStep(), e, myrank); + simulation.solve( GlobalSources[e], GlobalTimeSeries[e], mu, lambda, rho, U, Um, upred_saved, ucorr_saved, true, e, mopt->m_nsteps_in_memory, phcase, pseudo_hessian, step_to_record); sw4_profile->time_stamp("done forward solve" ); // Compute misfit, 'diffs' will hold the source for the adjoint problem @@ -609,7 +640,7 @@ void compute_f_and_df( EW& simulation, int nspar, int nmpars, double* xs, double dfsrc[11]; get_source_pars( nspar, dfsrc, dfs ); sw4_profile->time_stamp("backward+adjoint solve" ); - simulation.solve_backward_allpars( GlobalSources[e], rho, mu, lambda, diffs, U, Um, upred_saved, ucorr_saved, dfsrc, gRho, gMu, gLambda, e ); + simulation.solve_backward_allpars( GlobalSources[e], rho, mu, lambda, diffs, U, Um, upred_saved, ucorr_saved, dfsrc, gRho, gMu, gLambda, step_to_record, e ); sw4_profile->time_stamp("done backward+adjoint solve" ); // int ip=67, jp=20, kp=20; // if( simulation.interior_point_in_proc(ip,jp,0)) @@ -2123,11 +2154,11 @@ int main(int argc, char **argv) GlobalObservations, myRank, mopt ); else if( mopt->m_optmethod == 2 ) nlcg( simulation, nspar, nmpars, xs, nmpard, xm, GlobalSources, GlobalTimeSeries, - GlobalObservations, myRank, mopt ); - - for( int e=0 ; e < GlobalTimeSeries.size(); e++ ) + GlobalObservations, myRank, mopt ); + // save time windows + for( int e=0 ; e < GlobalTimeSeries.size(); e++ ) for( int m=0 ; m < GlobalTimeSeries[e].size() ; m++ ) - GlobalTimeSeries[e][m]->writeWindows(); + GlobalTimeSeries[e][m]->writeWindows(); sw4_profile->time_stamp("Done optimizer"); sw4_profile->flush(); diff --git a/src/solve-backward-allpars.C b/src/solve-backward-allpars.C index 078b127e..900b04d3 100644 --- a/src/solve-backward-allpars.C +++ b/src/solve-backward-allpars.C @@ -11,13 +11,14 @@ void EW::solve_backward_allpars( vector & a_Sources, vector & a_TimeSeries, vector& Up, vector& U, vector& Upred_saved, vector& Ucorr_saved, double gradientsrc[11], vector& gRho, vector& gMu, - vector& gLambda, int event ) + vector& gLambda, int step_to_record, int event ) { // solution arrays vector F, Lk, Kacc, Kp, Km, K, Um, Uacc; // vector gRho, gMu, gLambda; vector AlphaVE, AlphaVEm, AlphaVEp; vector BCForcing; + bool verbose=false; F.resize(mNumberOfGrids); Lk.resize(mNumberOfGrids); @@ -119,7 +120,10 @@ void EW::solve_backward_allpars( vector & a_Sources, K[g].set_to_zero(); } double t = mDt*(mNumberOfTimeSteps[event]-1); - int beginCycle = 1; + + //if(proc_zero() && verbose) cout << "first time step to record sides=" << step_to_record << endl; + + int beginCycle = 1; double time_measure[8]; double time_sum[8]={0,0,0,0,0,0,0,0}; @@ -198,7 +202,8 @@ void EW::solve_backward_allpars( vector & a_Sources, // set boundary data on Uacc, from forward solver for( int g=0 ; g < mNumberOfGrids ; g++ ) { - Upred_saved[g]->pop( Uacc[g], currentTimeStep ); + if(currentTimeStep>=step_to_record) Upred_saved[g]->pop( Uacc[g], currentTimeStep ); + else Uacc[g].set_to_zero(); communicate_array( Uacc[g], g ); } // Note, this assumes BCForcing is not time dependent, which is usually true @@ -210,7 +215,10 @@ void EW::solve_backward_allpars( vector & a_Sources, evalCorrector( Um, a_Rho, Lk, F ); for( int g=0 ; g < mNumberOfGrids ; g++ ) - Ucorr_saved[g]->pop( Um[g], currentTimeStep-2 ); + { + if(currentTimeStep>=step_to_record) Ucorr_saved[g]->pop( Um[g], currentTimeStep-2 ); + else Um[g].set_to_zero(); + } // set boundary data on U for(int g=0 ; g < mNumberOfGrids ; g++ ) diff --git a/src/solve.C b/src/solve.C index 0d72496c..2b2026f7 100644 --- a/src/solve.C +++ b/src/solve.C @@ -69,7 +69,7 @@ void EW::solve( vector & a_Sources, vector & a_TimeSeries, vector& U, vector& Um, vector& Upred_saved_sides, vector& Ucorr_saved_sides, bool save_sides, - int event, int nsteps_in_memory, int varcase, vector& PseudoHessian ) + int event, int nsteps_in_memory, int varcase, vector& PseudoHessian, int step_to_record) { // Experimental // int nsteps_in_memory=50; @@ -78,6 +78,7 @@ void EW::solve( vector & a_Sources, vector & a_TimeSeries, vector AlphaVE, AlphaVEm, AlphaVEp; // vectors of pointers to hold boundary forcing arrays in each grid vector BCForcing; + bool verbose=true; BCForcing.resize(mNumberOfGrids); F.resize(mNumberOfGrids); @@ -219,6 +220,9 @@ void EW::solve( vector & a_Sources, vector & a_TimeSeries, } } +// first step to record boundary values +// if(proc_zero() && verbose) cout << "solver: first time step to record sides=" << step_to_record << endl; + // Set the number of time steps, allocate the recording arrays, and set reference time in all time series objects /* #pragma omp parallel for */ for (int ts=0; ts & a_Sources, vector & a_TimeSeries, printf("End report of internal flags and settings\n\n"); } +/* +if(step_to_record==1) +{ if( save_sides ) { for( int g=0 ; g < mNumberOfGrids ; g++ ) { - Upred_saved_sides[g]->push( Um[g], -1 ); - Upred_saved_sides[g]->push( U[g], 0 ); - Ucorr_saved_sides[g]->push( Um[g], -1 ); - Ucorr_saved_sides[g]->push( U[g], 0 ); + Upred_saved_sides[g]->push( Um[g], -1 ); + Upred_saved_sides[g]->push( U[g], 0 ); + Ucorr_saved_sides[g]->push( Um[g], -1 ); + Ucorr_saved_sides[g]->push( U[g], 0 ); } } +} +*/ for( int g=0 ; g < mNumberOfGrids ; g++ ) Up[g].set_to_zero(); @@ -621,6 +630,17 @@ void EW::solve( vector & a_Sources, vector & a_TimeSeries, if( is_debug && proc_zero() ) cout << "Solving " << currentTimeStep << endl; +if( save_sides && currentTimeStep==(step_to_record>beginCycle? step_to_record : beginCycle)) + { + for( int g=0 ; g < mNumberOfGrids ; g++ ) + { + Upred_saved_sides[g]->push( Um[g], currentTimeStep-1-1 ); // save wavefields on the sides using push method of DataPatch + Upred_saved_sides[g]->push( U[g], currentTimeStep-1+0 ); + Ucorr_saved_sides[g]->push( Um[g], currentTimeStep-1-1 ); + Ucorr_saved_sides[g]->push( U[g], currentTimeStep-1+0 ); + } + } + // all types of forcing... bool trace =false; int dbgproc = 1; @@ -804,9 +824,11 @@ void EW::solve( vector & a_Sources, vector & a_TimeSeries, if( trace && m_myRank == dbgproc ) cout <<" after evalDpDmInTime" << endl; - if( save_sides ) - for( int g=0 ; g < mNumberOfGrids ; g++ ) - Upred_saved_sides[g]->push( Uacc[g], currentTimeStep ); + if( save_sides && currentTimeStep >= step_to_record) + { + for( int g=0 ; g < mNumberOfGrids ; g++ ) + Upred_saved_sides[g]->push( Uacc[g], currentTimeStep ); + } if( m_checkfornan ) check_for_nan( Uacc, 1, "uacc " ); @@ -926,9 +948,11 @@ void EW::solve( vector & a_Sources, vector & a_TimeSeries, if( m_do_geodynbc ) restore_geoghost(Up); - if( save_sides ) - for( int g=0 ; g < mNumberOfGrids ; g++ ) - Ucorr_saved_sides[g]->push( Up[g], currentTimeStep ); + if( save_sides && currentTimeStep >= step_to_record) + { + for( int g=0 ; g < mNumberOfGrids ; g++ ) + Ucorr_saved_sides[g]->push( Up[g], currentTimeStep ); + } }// end if mOrder == 4 diff --git a/src/solveTT.C b/src/solveTT.C index 87265004..74e8f766 100644 --- a/src/solveTT.C +++ b/src/solveTT.C @@ -10,6 +10,7 @@ #include "MaterialParameterization.h" #define STRINGSIZE 128 +#define FLOATMAX 1e20 #ifdef USE_HDF5 #include "sachdf5.h" @@ -17,6 +18,36 @@ #define SQR(x) ((x)*(x)) +// find minimum time at the surface +float time_to_record(float* t, int nx, int ny, int nz) +{ + float tmin=FLOATMAX; + unsigned long offset; + +// top and bottom + offset = (nz-1)*nx*ny; + for(int i=0; im_misfit == Mopt::CROSSCORR // output to a file for QC char file[STRINGSIZE]; - FILE *fd; + FILE *fd, *ft; if(myrank==0) { sprintf(file, "%s/eikonal_time_%d.txt", a_TimeSeries[0]->getPath().c_str(),event); fd = fopen(file, "w"); + sprintf(file, "%s/time_to_record_%d.txt", a_TimeSeries[0]->getPath().c_str(),event); + ft = fopen(file, "w"); + if(ft==NULL) fprintf(stderr, "unable to open file:%s\n", file); fprintf(fd, "event count station x y z tstart_vp tend_vp tstart_vs tend_vs\n"); } @@ -179,7 +213,10 @@ int ntr = a_TimeSeries.size(); } // loop over traces if(myrank==0) fclose(fd); - + if(myrank==0 && ft!=NULL) { + fprintf(ft, "%g\t%g\n", time_to_record(timep,nx,ny,nz), time_to_record(times,nx,ny,nz)); + fclose(ft); + } if(timep) free(timep); if(times) free(times); if(cp) free(cp); @@ -193,162 +230,5 @@ int ntr = a_TimeSeries.size(); } -/* -void EW::solveTT_old( Source* a_Source, vector & a_TimeSeries, - double* xs, int nmpars, MaterialParameterization* mp, int wave_mode, float win, float twinshift, int event, int myrank) -{ - int verbose=0; - int ix, iy, iz, nx, ny, nz, n; - nx = mp->getNX(); - ny = mp->getNY(); - nz = mp->getNZ(); - if( verbose >= 1 && myrank == 0 ) - std::cout << "source x0=" << a_Source->getX0() << " y0=" << a_Source->getY0() << " z0=" << a_Source->getZ0() << std::endl; - n= nmpars/nx/ny/nz; - if( verbose >= 1 && myrank == 0 ) - { - std::cout << "mp->xmin=" << mp->getXmin() << " ymin=" << mp->getYmin() << " zmin=" << mp->getZmin() << std::endl; - std::cout << "nx=" << nx << " ny=" << ny << " nz=" << nz << " n=" << n << std::endl; - } - -float *cp =(float*)malloc(nmpars/n*sizeof(float)); -float *cs =(float*)malloc(nmpars/n*sizeof(float)); - -float *sp =(float*)malloc(nmpars/n*sizeof(float)); -float *ss =(float*)malloc(nmpars/n*sizeof(float)); -float* timep =(float*)malloc(nmpars/n*sizeof(float)); -float* times =(float*)malloc(nmpars/n*sizeof(float)); -int *inflag =(int*)malloc(nmpars/n*sizeof(int)); - - -float_sw4 cpmin, cpmax; -float_sw4 csmin, csmax; - -cpmin=1e20; cpmax=-1e20; -csmin=1e20; csmax=-1e20; - - -//#pragma omp parallel for - for( size_t i = 0 ; i < nmpars/n ; i++ ) - { - cs[i] = xs[n*i+n-2]; - cp[i] = xs[n*i+n-1]; - - ss[i] = 1./xs[n*i+n-2]; - ss[i] = ss[i]*ss[i]; - sp[i] = 1./xs[n*i+n-1]; - sp[i] = sp[i]*sp[i]; - - if(cp[i]cpmax) cpmax=cp[i]; - if(cs[i]csmax) csmax=cs[i]; - } - - if( verbose >= 1 && myrank == 0 ) - { - std::cout << "solveTT cpmin=" << cpmin << " cpmax=" << cpmax << std::endl; - std::cout << "solveTT csmin=" << csmin << " csmax=" << csmax << std::endl; - } - -int ntr = a_TimeSeries.size(); -//std::cout << "nmpars=" << nmpars << " number of traces=" << ntr << std::endl; - - -//std::cout << "mp xmin=" << mp->getXmin() << " hx=" << mp->getDx() << " nx=" << mp->getNX() << std::endl; - bool plane[3]={false,false,false}; - - fastmarch_init(mp->getNZ(),mp->getNY(),mp->getNX()); // nx fast - fastmarch(timep, // time - sp, // slowness squared - inflag, // in/front/out flag - plane, // if plane source - nz, ny, nx, // dimensions - mp->getZmin(),mp->getYmin(),mp->getXmin(), // origin - mp->getDz(),mp->getDy(),mp->getDx(), // sampling - a_Source->getZ0(),a_Source->getY0(),a_Source->getX0(), // source - 1,1,1, // box around the source - 2); // accuracy order (1,2,3) - - fastmarch(times, // time - ss, // slowness squared - inflag, // in/front/out flag - plane, // if plane source - nz, ny, nx, // dimensions - mp->getZmin(),mp->getYmin(),mp->getXmin(), // origin - mp->getDz(),mp->getDy(),mp->getDx(), // sampling - a_Source->getZ0(),a_Source->getY0(),a_Source->getX0(), // source - 1,1,1, // box around the source - 2); // accuracy order (1,2,3) - - fastmarch_close(); - - //checkMinMax(nmpars/3, timep, "timep");mopt->m_misfit == Mopt::CROSSCORR - // output to a file for QC - char file[STRINGSIZE]; - FILE *fd; - - if(myrank==0) { - sprintf(file, "%s/time_event_%d.txt", a_TimeSeries[0]->getPath().c_str(),event); - fd = fopen(file, "w"); - } - - for(int ig=0; iggetX() - mp->getXmin())/mp->getDx()+0.5; - iy = (a_TimeSeries[ig]->getY() - mp->getYmin())/mp->getDy()+0.5; - iz = (a_TimeSeries[ig]->getZ() - mp->getZmin())/mp->getDz()+0.5; - - //std::cout << "trace ig=" << ig+1 << " x=" << a_TimeSeries[ig]->getX() << " y=" << a_TimeSeries[ig]->getY() << " tp=" - // << timep[iz*nx*ny+iy*nx+ix] << " ts=" << times[iz*nx*ny+iy*nx+ix] << std::endl; - - // set window - switch(wave_mode) { - case 0: // P-wave only - a_TimeSeries[ig]->set_window(timep[iz*nx*ny+iy*nx+ix]+win*twinshift, timep[iz*nx*ny+iy*nx+ix]+win, - timep[iz*nx*ny+iy*nx+ix]+win*twinshift, timep[iz*nx*ny+iy*nx+ix]+win); - - if(myrank==0) fprintf(fd, "%d %d\t%s\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", - event, ig+1, a_TimeSeries[ig]->getStationName().c_str(), a_TimeSeries[ig]->getX(),a_TimeSeries[ig]->getY(),a_TimeSeries[ig]->getZ(), - timep[iz*nx*ny+iy*nx+ix]+win*twinshift, timep[iz*nx*ny+iy*nx+ix]+win, - timep[iz*nx*ny+iy*nx+ix]+win*twinshift, timep[iz*nx*ny+iy*nx+ix]+win); - break; - case 1: // S-wave only - a_TimeSeries[ig]->set_window(times[iz*nx*ny+iy*nx+ix]+win*twinshift, times[iz*nx*ny+iy*nx+ix]+win, - times[iz*nx*ny+iy*nx+ix]+win*twinshift, times[iz*nx*ny+iy*nx+ix]+win); - - if(myrank==0) fprintf(fd, "%d %d\t%s\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", - event, ig+1, a_TimeSeries[ig]->getStationName().c_str(), a_TimeSeries[ig]->getX(),a_TimeSeries[ig]->getY(),a_TimeSeries[ig]->getZ(), - times[iz*nx*ny+iy*nx+ix]+win*twinshift, times[iz*nx*ny+iy*nx+ix]+win, - times[iz*nx*ny+iy*nx+ix]+win*twinshift, times[iz*nx*ny+iy*nx+ix]+win); - - break; - default: - a_TimeSeries[ig]->set_window(timep[iz*nx*ny+iy*nx+ix]+win*twinshift, timep[iz*nx*ny+iy*nx+ix]+win, - times[iz*nx*ny+iy*nx+ix]+win*twinshift, times[iz*nx*ny+iy*nx+ix]+win); - - if(myrank==0) fprintf(fd, "%d %d\t%s\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", - event, ig+1, a_TimeSeries[ig]->getStationName().c_str(), a_TimeSeries[ig]->getX(),a_TimeSeries[ig]->getY(),a_TimeSeries[ig]->getZ(), - timep[iz*nx*ny+iy*nx+ix]+win*twinshift, timep[iz*nx*ny+iy*nx+ix]+win, - times[iz*nx*ny+iy*nx+ix]+win*twinshift, times[iz*nx*ny+iy*nx+ix]+win); - } - } // loop over traces - - if(myrank==0) fclose(fd); - - - if(timep) free(timep); - if(times) free(times); - if(cp) free(cp); - if(cs) free(cs); - if(sp) free(sp); - if(ss) free(ss); - if(inflag) free(inflag); - - if( verbose >= 1 && myrank == 0 ) - std::cout << "solveTT completed" << std::endl; - -} -*/