// asg_scal_p.cc
// Author: Mario J. Bencomo
// last modified: 12/12/16
#include "asg_defn.hh"
#include "asg.hh"
#include "MPS_includes.hh"
#include "MPS_Space_Examples.hh"
#include "MPS_to_RHS.hh"
using RVL::valparse;
using RVL::AssignFilename;
using RVL::OperatorEvaluation;
using TSOpt::Scal_MPS_Space;
using TSOpt::MPS_to_RHS;
using TSOpt::add_to_pars;
const char *sdoc[] = {
" ============================================================================",
" asg_scal_p.x ",
" ============================================================================",
" Authors: Mario J. Bencomo",
" ",
" Acoustic variable density centered difference modeling for single scalar",
" multipole sources and pressure data output.",
" ",
" PDE:",
" dp/dt + bulk{ div v } = f(t,x)",
" dv/dt + buoy grad p = 0",
"",
" Source term:",
" 2-D, f(x,t) = f_0(t) (d/dx_0)^{s_0} (d/dx_1)^{s_1} delta(x-x^*)",
" 3-D, f(x,t) = f_0(t) (d/dx_0)^{s_0} (d/dx_1)^{s_1} (d/dx_2)^{s_2} delta(x-x^*)",
" ",
" Typical parameter list. May be copied, edited, and used for input: either",
" include parameters on command line (for example in Flow), or place",
" in file <foo> and include \"par=<foo>\" on command line. Any parameter",
" included explicitly in command line overrides parameter with same key",
" in par file.",
" ",
" Invoke single threaded execution by ",
" \"asg_p_scal.x [parameters] [standalone install]\"",
" ",
" or multi-threaded execution using interactive or batch MPI (for which",
" code must be compiled with MPI enabled).",
" ",
" Given values are defaults; non-optional values indicated by corner brackets.",
" ",
" --------------------------- begin parameters ---------------------------",
" FD info:",
" ",
" order = 2 spatial half-order",
" cfl = 0.75 proportion of max dt/dx",
" cmin = 1.0 min permitted velocity (m/ms) - sanity check",
" cmax = 4.5 max permitted velocity (m/ms) - used in dt comp",
" dmin = 1.0 min permitted density (g/cm^3) - sanity check",
" dmax = 1.0 max permitted density (g/cm^3) - sanity check",
" ",
" ------------------------------------------------------------------------",
" Source info:",
" ",
" source_p=empty dummy key required for asg driver.",
" ",
" ------------------------------------------------------------------------",
" Trace info:",
" ",
" data_p = <path> data file, SU format - must exist ",
" headers establish acquisition geometry",
" output for forward modeling (including all",
" orders of derivative), input for adjoint",
" ",
" ------------------------------------------------------------------------",
" Model info:",
" ",
" bulkmod = <path> input bulk modulus, ",
" determines simulation spatial grid,",
" buoyancy= <path> input buoyancy",
" ",
" ------------------------------------------------------------------------",
" MPI info:",
" ",
" mpi_np1 = 1 number of subdomains along axis 1",
" mpi_np2 = 1 number of subdomains along axis 2",
" mpi_np3 = 1 number of subdomains along axis 3",
" partask = 1 number of shots to execute in parallel",
" ",
" ------------------------------------------------------------------------",
" Output info:",
" ",
" FD ouput - written to coutxxx.txt on rank xxx",
" printact = 0 monitor output:",
" < 0 - none",
" 0 - announce each shot simulation",
" 1 - print time step index",
" 2 - diagnostic messages from main routines",
" > 2 - much more, useful only for serious ",
" debugging",
" dump_pi = 0 dump parallel/dom. decomp info",
" dump_lda = 0 dump grid data for allocated arrays",
" dump_ldc = 0 dump grid data for computational arrays",
" dump_ldr = 0 dump grid data for receive buffers",
" dump_lds = 0 dump grid data for send buffers",
" dump_term = 0 dump trace header data",
" ",
" ------------------------------------------------------------------------",
" MPS info:",
" ",
" appx_ord = 4 singular source approximation order ",
" should match spatial order of FD method, i.e., 2*order",
" MPS_ord_0 = 0 MPS order in z-direction, s_0",
" MPS_ord_1 = 0 MPS order in x-direction, s_1",
" MPS_ord_2 = 0 MPS order in y-direction, s_2",
" MPS_file = <path> MPS file containing MPS coefficients, SU format - must exist",
" Must contain source locations.",
" ---------------------------end parameters ------------------------------",
NULL };
IOKEY IWaveInfo::iwave_iokeys[]
= {
{"bulkmod", 0, true, 1 },
{"buoyancy", 1, true, 1 },
{"source_p", 2, true, 2 },
{"data_p", 2, false, 2 },
{"", 0, false, 0 }
};
int xargc;
char **xargv;
int main(int argc, char ** argv) {
try {
#ifdef IWAVE_USE_MPI
int ts=0;
MPI_Init_thread(&argc,&argv,MPI_THREAD_FUNNELED,&ts);
#endif
PARARRAY * pars = NULL;
FILE * stream = NULL;
TSOpt::IWaveEnvironment(argc, argv, 0, &pars, &stream);
#ifdef IWAVE_USE_MPI
if(retrieveGlobalRank()==0){
#endif
if(argc<2){
pagedoc();
exit(0);
}
cerr << "\n////////////////////////////////////////////////////////////\n"
<< "Running asg_scal_p.x\n";
#ifdef IWAVE_USE_MPI
}
#endif
TSOpt::MPS_KEYS mks;
mks.appx_ord = "appx_ord";
mks.MPS_file = "MPS_file";
mks.grid_file = "bulkmod";
mks.RHS_files.push_back("source_p");
//initializing MPS space
Scal_MPS_Space mps_sp(mks,*pars);
//initializing MPS_to_RHS
MPS_to_RHS m2r(mks,*pars,mps_sp);
//constructing MPS IWaveLOVOp
TSOpt::IWaveLOVOp lovop(m2r.get_pars(),stream);
RVL::LinCompLOVOp<float> mps_lovop(m2r,lovop);
//setting up domain vector
Vector<ireal> x(mps_lovop.getDomain());
Components<ireal> cx(x);
Components<ireal> cx_m(cx[0]); //model params
RVL::AssignFilename af_bmd(valparse<string>(*pars,"bulkmod"));
cx_m[0].eval(af_bmd);
RVL::AssignFilename af_buo(valparse<string>(*pars,"buoyancy"));
cx_m[1].eval(af_buo);
Components<ireal> cx_f(cx[1]); //source params
RVL::AssignFilename af_mps(valparse<string>(*pars,"MPS_file"));
cx_f[0].eval(af_mps);
//setting up range vector
Vector<ireal> y(mps_lovop.getRange());
RVL::AssignFilename af_data(RVL::valparse<string>(*pars,"data_p"));
y.eval(af_data);
//applying op
RVL::OperatorEvaluation<ireal> eval(mps_lovop,x);
y.copy(eval.getValue());
//clean up
ps_delete(&pars);
iwave_fdestroy();
#ifdef IWAVE_USE_MPI
if(retrieveGlobalRank()==0){
#endif
cerr << "\nFinishing asg_scal_p.x\n"
<< "////////////////////////////////////////////////////////////\n";
#ifdef IWAVE_USE_MPI
}
MPI_Finalize();
#endif
}
catch(RVLException &e){
e << "Exiting with error!\n";
#ifdef IWAVE_USE_MPI
if( retrieveGlobalRank()==0 ){
#endif
e.write(cerr);
#ifdef IWAVE_USE_MPI
}
MPI_Barrier(retrieveGlobalComm());
MPI_Abort(retrieveGlobalComm(),0);
#endif
exit(1);
}
}