Skip to content

Memory/computation enhancement through bypassing recording of sides over source precursor #129

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 8 commits into
base: developer
Choose a base branch
from
Draft

Memory/computation enhancement through bypassing recording of sides over source precursor #129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
0d430df
draft on skipping t0 in inversion
Feb 4, 2022
3253980
remove time shifting and add t0 to time window
Feb 16, 2022
c5d0435
Merge branch 'wl/addt0' into wl/skipoffset
Feb 23, 2022
33732a4
add input option to turn on/off skipping recording sides over source …
Feb 25, 2022
d1575b5
Merge branch 'developer' into wl/skipoffset
Mar 1, 2022
c658f95
correct first instance of recording at sides
Mar 1, 2022
248ef9e
revise backward solver to reach the end, add an estimate of minimum t…
Mar 22, 2022
7549dd6
minor print correction
Mar 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/EW.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ void solve( vector<Source*> & a_GlobalSources, vector<TimeSeries*> & a_GlobalTim
vector<Sarray>& U, vector<Sarray>& Um,
vector<DataPatches*>& Upred_saved_sides,
vector<DataPatches*>& Ucorr_saved_sides, bool save_sides, int event, int save_steps,
int varcase, vector<Sarray>& pseudoHessian );
int varcase, vector<Sarray>& pseudoHessian, int step_to_record);

void solveTT( Source* a_GlobalSource, vector<TimeSeries*> & a_GlobalTimeSeries,
double* xs, int nmpars, MaterialParameterization* mp, int wave_mode, int event, int myrank);
Expand All @@ -134,7 +134,7 @@ void solve_backward_allpars( vector<Source*> & a_GlobalSources, vector<Sarray>&
vector<Sarray>& a_Lambda, vector<TimeSeries*> & a_GlobalTimeSeries,
vector<Sarray>& a_U, vector<Sarray>& a_Um, vector<DataPatches*>& Upred_saved_sides,
vector<DataPatches*>& Ucorr_saved_sides, float_sw4 gradients[11],
vector<Sarray>& gRho, vector<Sarray>& gMu, vector<Sarray>& gLambda, int event );
vector<Sarray>& gRho, vector<Sarray>& gMu, vector<Sarray>& gLambda, int step_to_record, int event );
//int nmpar, float_sw4* gradientm );

void solve_dudp( vector<Source*>& a_Sources, vector<Sarray>& a_Rho,
Expand Down
8 changes: 8 additions & 0 deletions src/Mopt.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
2 changes: 2 additions & 0 deletions src/Mopt.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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; }
Expand All @@ -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;
Expand Down
6 changes: 3 additions & 3 deletions src/Source.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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();

Expand Down Expand Up @@ -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();

Expand Down
1 change: 0 additions & 1 deletion src/TimeSeries.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
12 changes: 11 additions & 1 deletion src/lbfgs.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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() )
Expand All @@ -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 )
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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<maxit && !converged


Expand Down
2 changes: 1 addition & 1 deletion src/main.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ main(int argc, char **argv)
vector<Sarray> 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

Expand Down
45 changes: 38 additions & 7 deletions src/moptmain.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -225,6 +225,33 @@ void set_timewindows_from_eikonal_time(vector<vector<TimeSeries*> >& GlobalTimeS
} // end of events
}

int firstTimeStep_to_recordBC(const vector<Source*> & a_Sources, const vector<TimeSeries*> & 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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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();
Expand Down
16 changes: 12 additions & 4 deletions src/solve-backward-allpars.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,14 @@ void EW::solve_backward_allpars( vector<Source*> & a_Sources,
vector<TimeSeries*> & a_TimeSeries, vector<Sarray>& Up, vector<Sarray>& U,
vector<DataPatches*>& Upred_saved, vector<DataPatches*>& Ucorr_saved,
double gradientsrc[11], vector<Sarray>& gRho, vector<Sarray>& gMu,
vector<Sarray>& gLambda, int event )
vector<Sarray>& gLambda, int step_to_record, int event )
{
// solution arrays
vector<Sarray> F, Lk, Kacc, Kp, Km, K, Um, Uacc;
// vector<Sarray> gRho, gMu, gLambda;
vector<Sarray*> AlphaVE, AlphaVEm, AlphaVEp;
vector<double **> BCForcing;
bool verbose=false;

F.resize(mNumberOfGrids);
Lk.resize(mNumberOfGrids);
Expand Down Expand Up @@ -119,7 +120,10 @@ void EW::solve_backward_allpars( vector<Source*> & 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};
Expand Down Expand Up @@ -198,7 +202,8 @@ void EW::solve_backward_allpars( vector<Source*> & 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
Expand All @@ -210,7 +215,10 @@ void EW::solve_backward_allpars( vector<Source*> & 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++ )
Expand Down
46 changes: 35 additions & 11 deletions src/solve.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & a_TimeSeries,
vector<Sarray>& U, vector<Sarray>& Um,
vector<DataPatches*>& Upred_saved_sides,
vector<DataPatches*>& Ucorr_saved_sides, bool save_sides,
int event, int nsteps_in_memory, int varcase, vector<Sarray>& PseudoHessian )
int event, int nsteps_in_memory, int varcase, vector<Sarray>& PseudoHessian, int step_to_record)
{
// Experimental
// int nsteps_in_memory=50;
Expand All @@ -78,6 +78,7 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & a_TimeSeries,
vector<Sarray*> AlphaVE, AlphaVEm, AlphaVEp;
// vectors of pointers to hold boundary forcing arrays in each grid
vector<float_sw4**> BCForcing;
bool verbose=true;

BCForcing.resize(mNumberOfGrids);
F.resize(mNumberOfGrids);
Expand Down Expand Up @@ -219,6 +220,9 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & 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_TimeSeries.size(); ts++)
Expand Down Expand Up @@ -579,16 +583,21 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & 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();
Expand Down Expand Up @@ -621,6 +630,17 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & 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;
Expand Down Expand Up @@ -804,9 +824,11 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & 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 " );
Expand Down Expand Up @@ -926,9 +948,11 @@ void EW::solve( vector<Source*> & a_Sources, vector<TimeSeries*> & 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

Expand Down
Loading