-
JanThorbecke authoredJanThorbecke authored
getParameters.c 48.88 KiB
#include<stdlib.h>
#include<stdio.h>
#include<stdbool.h>
#include<math.h>
#include<float.h>
#include<time.h>
#include"fdacrtmc.h"
#include"par.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
int readModel(modPar *mod);
int readDT(srcPar *src, modPar *mod);
int createRcvCoordinates(modPar *mod,srcPar *rcv, recPar *rec, int verbose);
int prepareFDOperators(modPar *mod, bndPar *bnd, decompPar *decomp);
int setWisdomPath(char *path);
void printWisdomPath();
int initFFTW(void);
int initPML(modPar *mod, bndPar *bnd, int verbose);
int readRcvCoordinates(modPar *mod, srcPar *rcv, recPar *rec, int verbose);
int initializeCompressionField(migPar *mig);
int getParameters(modPar *mod, srcPar *src, srcPar *rcv, bndPar *bnd, snaPar *sna, migPar *mig, recPar *rec, decompPar *decomp, int *verbose){
float cmax, cmin, scl, wfct, dispfactor, stabfactor;
int tmpint;
size_t ix, iz, ioPz;
/***************************/
/* Input Parameters Part 1 */
/***************************/
// Verbosity
if(!getparint("verbose",verbose)) *verbose=0; //Default: Silent
// Modelling Scheme
if(!getparint("ischeme",&mod->ischeme)) mod->ischeme=1; //Default: Acoustic
if(!getparint("iorder",&mod->iorder)) mod->iorder=4; //Default: 4th Order
if(!getparint("sh",&mod->sh)) mod->sh=0; //???? SH-Waves
// Model Files
if(!getparstring("file_cp",&mod->file_cp)) verr("P-Wave Velocity File Required!");
if(mod->ischeme>2) if(!getparstring("file_cs",&mod->file_cs)) verr("S-Wave Velocity File Required For Elastic Modeling!");
if(!getparstring("file_den",&mod->file_den)) mod->file_den=NULL;
// Input Wavefield Files
if(!getparstring("file_src" ,&src->file_src)) verr("Source File For Forward Propagation Required!");
if(!getparstring("file_rcv" ,&rcv->file_src)) rcv->file_src=NULL;
if(!getparstring("file_loc" ,&rec->file_loc)) rec->file_loc=NULL;
if(!getparint( "rcv_left" ,&rec->left )) rec->left =0;
if(!getparint( "rcv_top" ,&rec->top )) rec->top =0;
if(!getparint( "rcv_bottom",&rec->bottom )) rec->bottom=0;
if(!getparint( "rcv_right" ,&rec->right )) rec->right =0;
if(!getparint( "rcv_vx" ,&rec->vx )) rec->vx =0;
if(!getparint( "rcv_vz" ,&rec->vz )) rec->vz =0;
if(mod->ischeme<2){
if(!getparint("rcv_p" ,&rec->p)) rec->p =0;
rec->txx=0;rec->txz=0;rec->tzz=0;
}else{
rec->p=0;
if(!getparint("rcv_txx" ,&rec->txx)) rec->txx=0;
if(!getparint("rcv_txz" ,&rec->txz)) rec->txz=0;
if(!getparint("rcv_tzz" ,&rec->tzz)) rec->tzz=0;
}
rec->rec=0; if(rec->left||rec->top||rec->right||rec->bottom){
if(mod->ischeme<2||rec->p) rec->rec=1; //Acoustic
else if(rec->txx||rec->txz||rec->tzz) rec->rec=1; //Elastic
if(rec->vx||rec->vz) rec->rec=1;
}
if(!getparint("rcv_write" ,&rec->write)) rec->write =0;
if(rec->write&&!rcv->file_src){
vwarn("No Receiver File Defined For Writing!");
rec->write=0;
}
// FFTw Wisdom Directory
if(getparstring("file_WisdomPath",&mig->file_mig)) setWisdomPath(mig->file_mig);
// Initialize Multithreaded FFTw
tmpint=initFFTW();
if(tmpint==-1) verr("Multithreaded FFTw Wisdom Stem Too Long");
else if(tmpint) verr("Failed to Initialize Multithreaded FFTw");
// Output Files
if(!getparstring("file_mig",&mig->file_mig)) mig->file_mig="mig.su";
if(!getparstring("file_snap",&sna->file_snap)) sna->file_snap="snap.su";
if(!getparstring("file_imp",&mod->file_imp)) mod->file_imp=NULL;
if(!getparstring("file_dd",&mod->file_dd)) mod->file_dd=NULL;
if(!getparint("writeDD",&decomp->writeDD)) decomp->writeDD=0;
if(decomp->writeDD&&!(*mod->file_dd)){
vwarn("Specified to write decomposition direction to disk, but no file name given.");
vwarn("Not writing decomposition direction to disk.");
decomp->writeDD=0;
}
/***************/
/* Read Models */
/***************/
readModel(mod);
/********************************/
/* Compute Statistics On Models */
/********************************/
// Max/Min Density
cmax=0;cmin=FLT_MAX;
#pragma omp parallel for reduction(max:cmax) reduction(min:cmin) private (ix,iz)
for(ix=0;ix<mod->nx;ix++){
#pragma ivdep
for(iz=0;iz<mod->nz;iz++){
if(cmax<mod->rho[ix*mod->nz+iz]) cmax=mod->rho[ix*mod->nz+iz];
if(cmin>mod->rho[ix*mod->nz+iz]) cmin=mod->rho[ix*mod->nz+iz];
}
}
mod->rho_max=cmax;mod->rho_min=cmin;
// Max/Min P-Wave Velocity
cmax=0;cmin=FLT_MAX;
#pragma omp parallel for reduction(max:cmax) reduction(min:cmin) private (ix,iz)
for(ix=0;ix<mod->nx;ix++){
#pragma ivdep
for(iz=0;iz<mod->nz;iz++){
if(cmax<mod->cp[ix*mod->nz+iz]) cmax=mod->cp[ix*mod->nz+iz];
if(cmin>mod->cp[ix*mod->nz+iz]) cmin=mod->cp[ix*mod->nz+iz];
}
}
mod->cp_max=cmax;mod->cp_min=cmin;
if(mod->ischeme>2){ /* Elastic */
// Max/Min P-Wave Velocity
cmax=0;cmin=FLT_MAX;
#pragma omp parallel for reduction(max:cmax) reduction(min:cmin) private (ix,iz)
for(ix=0;ix<mod->nx;ix++){
#pragma ivdep
for(iz=0;iz<mod->nz;iz++){
if(cmax<mod->cs[ix*mod->nz+iz]) cmax=mod->cs[ix*mod->nz+iz];
if(cmin>mod->cs[ix*mod->nz+iz]) cmin=mod->cs[ix*mod->nz+iz]; }
}
mod->cs_max=cmax;mod->cs_min=cmin;
cmax=MAX(mod->cp_max,mod->cs_max);
cmin=MIN(mod->cp_min,mod->cs_min);
}
/***************************/
/* Input Parameters Part 2 */
/***************************/
// Boundary Information
/* 1=free 2=pml 3=rigid 4=taper */
if(!getparint("left",&bnd->lef))bnd->lef=2;
if(!getparint("right",&bnd->rig))bnd->rig=2;
if(!getparint("top",&bnd->top))bnd->top=1;
if(!getparint("bottom",&bnd->bot))bnd->bot=2;
bnd->ntapo=0;
if(bnd->lef==2||bnd->rig==2||bnd->top==2||bnd->bot==2){
if(!getparint("npml",&sna->beam)) verr("For PML Boundaries Boundary Size Is Required");
bnd->npml=(size_t)sna->beam;
}else bnd->npml=0;
if(bnd->lef==4||bnd->rig==4||bnd->top==4||bnd->bot==4){
if(!getparint("ntaper",&sna->beam)) verr("For Tapered Boundaries Boundary Size Is Required");
bnd->ntapo=(size_t)sna->beam;
}else bnd->ntapo=0;
if(!getparfloat("R",&bnd->R))bnd->R=1e-5;
if(!getparfloat("m",&bnd->m))bnd->m=2.0;
if(bnd->ntapo<bnd->npml){bnd->ntap=bnd->npml;}else{bnd->ntap=bnd->ntapo;}
if(bnd->ntap!=0) if(!getparfloat("tapfact",&bnd->tapfact)) bnd->tapfact=0.30;
/********************************/
/* Process Boundary Information */
/********************************/
if(bnd->ntap){
bnd->tapx =(float*)malloc(bnd->ntap*sizeof(float));
bnd->tapz =(float*)malloc(bnd->ntap*sizeof(float));
bnd->tapxz=(float*)malloc(bnd->ntap*bnd->ntap*sizeof(float));
scl=bnd->tapfact/((float)bnd->ntap);
for(ix=0;ix<bnd->ntap;ix++){
wfct=(scl*ix);
bnd->tapx[ix]=exp(-(wfct*wfct));
wfct=(scl*(ix+0.5));
bnd->tapz[ix]=exp(-(wfct*wfct));
}
for(ix=0;ix<bnd->ntap;ix++){
for(iz=0;iz<bnd->ntap;iz++) {
wfct=(scl*sqrt(iz*iz+ix*ix));
bnd->tapxz[ix*bnd->ntap+iz]=exp(-(wfct*wfct));
}
}
}
/* Total number of gridpoints */
mod->naz=mod->nz+mod->iorder/2+1;
mod->nax=mod->nx+mod->iorder/2+1;
/* Vx: rox */
mod->ioXxb=mod->iorder/2;
mod->ieXxb=mod->nx+mod->ioXxb-1;
mod->ioXzb=mod->iorder/2-1;
mod->ieXzb=mod->nz+mod->ioXzb;
/* Vz: roz */
mod->ioZxb=mod->iorder/2-1;
mod->ieZxb=mod->nx+mod->ioZxb;
mod->ioZzb=mod->iorder/2;
mod->ieZzb=mod->nz+mod->ioZzb-1;
/* P, Txx, Tzz: lam, l2m */
mod->ioPxb=mod->iorder/2-1;
mod->iePxb=mod->nx+mod->ioPxb;
mod->ioPzb=mod->ioPxb;
mod->iePzb=mod->nz+mod->ioPzb; /* Txz: mul */
mod->ioTxb=mod->iorder/2;
mod->ieTxb=mod->nx+mod->ioTxb;
mod->ioTzb=mod->ioTxb;
mod->ieTzb=mod->nz+mod->ioTzb;
/* Add additional points for PML & tapered boundaries */
if(bnd->top==2||bnd->top==4){ // Top Boundary
mod->naz +=bnd->ntap;
// Horizontal Particle Velocity Grid
mod->ioXz=mod->ioXzb+bnd->ntap;
mod->ieXzb+=bnd->ntap;
// Vertical Particle Velocity Grid
mod->ioZz=mod->ioZzb+bnd->ntap;
mod->ieZzb+=bnd->ntap;
// Pressure Grid
mod->ioPz=mod->ioPzb+bnd->ntap;
mod->iePzb+=bnd->ntap;
// Txz Grid
mod->ioTz=mod->ioTzb+bnd->ntap;
mod->ieTzb+=bnd->ntap;
}else{
// Horizontal Particle Velocity Grid
mod->ioXz=mod->ioXzb;
// Vertical Particle Velocity Grid
mod->ioZz=mod->ioZzb;
// Pressure Grid
mod->ioPz=mod->ioPzb;
// Txz Grid
mod->ioTz=mod->ioTzb;
}
mod->ieXz=mod->ieXzb;
mod->ieZz=mod->ieZzb;
mod->iePz=mod->iePzb;
mod->ieTz=mod->ieTzb;
if(bnd->bot==2||bnd->bot==4){ // Bottom Boundary
mod->naz +=bnd->ntap;
// Horizontal Particle Velocity Grid
mod->ieXzb+=bnd->ntap;
// Vertical Particle Velocity Grid
mod->ieZzb+=bnd->ntap;
// Pressure Grid
mod->iePzb+=bnd->ntap;
// Txz Grid
mod->ieTzb+=bnd->ntap;
}
if(bnd->lef==2||bnd->lef==4){ // Left Boundary
mod->nax +=bnd->ntap;
// Horizontal Particle Velocity Grid
mod->ioXx=mod->ioXxb+bnd->ntap;
mod->ieXxb+=bnd->ntap;
// Vertical Particle Velocity Grid
mod->ioZx=mod->ioZxb+bnd->ntap;
mod->ieZxb+=bnd->ntap;
// Pressure Grid
mod->ioPx=mod->ioPxb+bnd->ntap;
mod->iePxb+=bnd->ntap;
// Txz Grid
mod->ioTx=mod->ioTxb+bnd->ntap;
mod->ieTxb+=bnd->ntap;
}else{
// Horizontal Particle Velocity Grid
mod->ioXx=mod->ioXxb;
// Vertical Particle Velocity Grid
mod->ioZx=mod->ioZxb;
// Pressure Grid
mod->ioPx=mod->ioPxb;
// Txz Grid
mod->ioTx=mod->ioTxb;
} mod->ieXx=mod->ieXxb;
mod->ieZx=mod->ieZxb;
mod->iePx=mod->iePxb;
mod->ieTx=mod->ieTxb;
if(bnd->rig==2||bnd->rig==4){ // Right Boundary
mod->nax +=bnd->ntap;
// Horizontal Particle Velocity Grid
mod->ieXxb+=bnd->ntap;
// Vertical Particle Velocity Grid
mod->ieZxb+=bnd->ntap;
// Pressure Grid
mod->iePxb+=bnd->ntap;
// Txz Grid
mod->ieTxb+=bnd->ntap;
}
mod->sizem=mod->nax*mod->naz; //Full Model Size
/* Initialize the array which contains the topography surface */
if(bnd->top==2||bnd->top==4){ioPz=mod->ioPzb;}else{ioPz=mod->ioPz;}
bnd->surface=(size_t*)malloc((mod->nax+mod->naz)*sizeof(size_t));
for(ix=0;ix<mod->nax+mod->naz;ix++)bnd->surface[ix]=ioPz;
if((bnd->top==2)||(bnd->bot==2)||(bnd->lef==2)||(bnd->rig==2)) bnd->pml=1;else bnd->pml=0;
/*********************************/
/* Initialize Receiver Locations */
/*********************************/
rcv->nsrc=0;
// On Boundary
if(rec->rec)createRcvCoordinates(mod,rcv,rec,*verbose);
// From File
if(rec->file_loc){readRcvCoordinates(mod,rcv,rec,*verbose);rec->rec=1;}
/*******************************/
/* Read Temporal Sampling Rate */
/*******************************/
readDT(src,mod);
/***************************/
/* Input Parameters Part 3 */
/***************************/
// Snapshot Information
if(!getparfloat("tsnap1" ,&sna->tsnap1 ))sna->tsnap1=0.1;
if(!getparfloat("tsnap2" ,&sna->tsnap2 ))sna->tsnap2=0.0;
if(!getparfloat("dtsnap" ,&sna->dt ))sna->dt=0.1 ;
if(!getparfloat("dxsnap" ,&sna->dx ))sna->dx=mod->dx;
if(!getparfloat("dzsnap" ,&sna->dz ))sna->dz=mod->dz;
if(!getparint( "snapwithbnd",&sna->withbnd))sna->withbnd=0 ;
if(sna->withbnd){
if(!getparfloat("xsnap1",&sna->xsnap1))sna->xsnap1=mod->origx-(mod->ioPx+1)*mod->dx ;
if(!getparfloat("xsnap2",&sna->xsnap2))sna->xsnap2=mod->xmax+(mod->nax-1-mod->iePx)*mod->dx;
if(!getparfloat("zsnap1",&sna->zsnap1))sna->zsnap1=mod->origz-(mod->ioPz+1)*mod->dz ;
if(!getparfloat("zsnap2",&sna->zsnap2))sna->zsnap2=mod->zmax+(mod->naz-1-mod->iePz)*mod->dz;
}else{
if(!getparfloat("xsnap1",&sna->xsnap1))sna->xsnap1=mod->origx;
if(!getparfloat("xsnap2",&sna->xsnap2))sna->xsnap2=mod->xmax ;
if(!getparfloat("zsnap1",&sna->zsnap1))sna->zsnap1=mod->origz;
if(!getparfloat("zsnap2",&sna->zsnap2))sna->zsnap2=mod->zmax ;
}
if(!getparint("snap_forw" ,&sna->forw ))sna->forw =1;
if(!getparint("snap_back" ,&sna->back ))sna->back =1;
if(!getparint("snap_vxshift" ,&sna->vxshift ))sna->vxshift=0;
if(!getparint("snap_vzshift" ,&sna->vzshift ))sna->vzshift=0;
if(!getparint("snap_vxtime" ,&sna->vxtime ))sna->vxtime =0;
if(!getparint("snap_vztime" ,&sna->vztime ))sna->vztime =0;
if(!getparint("beam" ,&sna->beam ))sna->beam =0;
if(!getparint("snap_type_vx" ,&sna->type.vx ))sna->type.vx=0;
if(!getparint("snap_type_vz" ,&sna->type.vz ))sna->type.vz=0;
if(mod->ischeme<3){ if(!getparint("snap_type_pu",&sna->type.pu))sna->type.pu=0;
if(!getparint("snap_type_pd",&sna->type.pd))sna->type.pd=0;
if(!getparint("snap_type_pl",&sna->type.pl))sna->type.pl=0;
if(!getparint("snap_type_pr",&sna->type.pr))sna->type.pr=0;
if(!getparint("snap_type_pn",&sna->type.pn))sna->type.pn=0;
sna->type.p=0;
if(!getparint("snap_type_p" ,&sna->type.p )&&!sna->type.vx&&!sna->type.vz&&!sna->type.pu&&!sna->type.pd&&!sna->type.pl&&!sna->type.pr&&!sna->type.pn)sna->type.p=1;
if(!sna->type.vx&&!sna->type.vz&&!sna->type.p&&!sna->type.pu&&!sna->type.pd&&!sna->type.pl&&!sna->type.pr&&!sna->type.pn){
vwarn("tsnap2>tsnap1 but no snapshot type selected. Not writing out snapshots.");
sna->tsnap2=0.0;
}
sna->type.txx=0;sna->type.txz=0;sna->type.tzz=0;
sna->type.pp =0;sna->type.ss =0;
}else if(mod->ischeme<5){
if(!getparint("snap_type_txx",&sna->type.txx))sna->type.txx=0;
if(!getparint("snap_type_txz",&sna->type.txz))sna->type.txz=0;
if(!getparint("snap_type_tzz",&sna->type.tzz))sna->type.tzz=0;
if(!getparint("snap_type_pp" ,&sna->type.pp ))sna->type.pp =0;
if(!getparint("snap_type_ss" ,&sna->type.ss ))sna->type.ss =0;
if(!sna->type.vx&&!sna->type.vz&&!sna->type.txx&&!sna->type.txz&&!sna->type.tzz&&!sna->type.pp&&!sna->type.ss){
vwarn("tsnap2>tsnap1 but no snapshot type selected. Not writing out snapshots.");
sna->tsnap2=0.0;
}
sna->type.p =0;sna->type.pu=0;sna->type.pd=0;sna->type.pl=0;
sna->type.pr=0;sna->type.pn=0;
}
// Migration Image Information
if(!getparint( "mig_mode", &mig->mode ))mig->mode=1 ;
if(!getparint( "mig_orient", &mig->orient ))mig->orient=1 ;
if(!getparint( "mig_backscatter",&mig->backscatter))mig->backscatter=0 ;
if(!getparint( "mig_tap", &mig->tap ))mig->tap=0 ;
if(!getparfloat("migdt", &mig->dt ))mig->dt=mod->dt ;
if(!getparfloat("migdx", &mig->dx ))mig->dx=mod->dx ;
if(!getparfloat("migdz", &mig->dz ))mig->dz=mod->dz ;
if(!getparfloat("migx1", &mig->xmig1 ))mig->xmig1=mod->origx;
if(!getparfloat("migx2", &mig->xmig2 ))mig->xmig2=mod->xmax ;
if(!getparfloat("migz1", &mig->zmig1 ))mig->zmig1=mod->origz;
if(!getparfloat("migz2", &mig->zmig2 ))mig->zmig2=mod->zmax ;
//TODO: Test Below
if(!getparint("mig_pp",&tmpint))tmpint=0;if(tmpint)mig->pp=true;else mig->pp=false;
if(!getparint("mig_pm",&tmpint))tmpint=0;if(tmpint)mig->pm=true;else mig->pm=false;
if(!getparint("mig_mp",&tmpint))tmpint=0;if(tmpint)mig->mp=true;else mig->mp=false;
if(!getparint("mig_mm",&tmpint))tmpint=0;if(tmpint)mig->mm=true;else mig->mm=false;
// Decomposition Information
decomp->direct=1;decomp->med=1;
if(!getparint("decomp_mavgi",&mod->mavgi ))mod->mavgi=0;
if(!getparint("decomp_mavga",&decomp->mavga )) decomp->mavga=0;
if(!getparint("decomp_mavgn",&decomp->mavgn )){decomp->mavgn=5,decomp->direct=0;}
if(!getparint("decomp_mavgo",&decomp->mavgo )){decomp->mavgo=3,decomp->med =0;}
if( getparint("decomp_mavg" ,&decomp->wavFilt)){
if(decomp->wavFilt!=0&&decomp->wavFilt!=1&&decomp->wavFilt!=3&&decomp->wavFilt!=5&&decomp->wavFilt!=7&&decomp->wavFilt!=9){
vwarn("decomp_mavg=%d is an unknown option, ignoring! Valid options are 0,1,3,5,7,9.");
decomp->mavgn=5;decomp->mavgo=3;
}else{
if(decomp->direct){
vwarn("decomp_mavg=%d, but decomp_mavgn=%d already specified, ignoring decomp_mavg for decomp_mavgn!",decomp->wavFilt,decomp->mavgn);
if(decomp->med){
vwarn("decomp_mavg=%d, but decomp_mavgo=%d already specified, ignoring decomp_mavg for decomp_mavgo!",decomp->wavFilt,decomp->mavgo);
}else{
decomp->mavgo=decomp->wavFilt;
}
}else{
decomp->mavgn=decomp->wavFilt;
if(decomp->med){
vwarn("decomp_mavg=%d, but decomp_mavgo=%d already specified, ignoring decomp_mavg for decomp_mavgo!",decomp->wavFilt,decomp->mavgo);
}else{
decomp->mavgo=decomp->wavFilt-2;
if(decomp->mavgo<0)decomp->mavgo=0;
}
} }
}
if(!getparint( "decomp_direct" ,&decomp->direct ))decomp->direct =0 ;
if(!getparint( "decomp_med" ,&decomp->med ))decomp->med =3 ;
if(!getparint( "decomp_wavFilt",&decomp->wavFilt))decomp->wavFilt=0 ;
if(!getparfloat("decomp_reg" ,&decomp->reg ))decomp->reg =0.001;
if(!getparfloat("decomp_kl" ,&decomp->kl ))decomp->kl =1.0 ;
if(!getparfloat("decomp_kh" ,&decomp->kh ))decomp->kh =1.5 ;
if(!getparfloat("decomp_px" ,&decomp->px ))decomp->px =0.0 ;
if(!getparfloat("decomp_pz" ,&decomp->pz ))decomp->pz =0.0 ;
if(!decomp->px&&!decomp->pz)decomp->pz=1.0;
if(decomp->direct){decomp->decomp=1,decomp->wavFilt=1;}
// Additional Information
if(mod->ischeme>2){
sna->type.p=0;
if(!getparint("sna_type_txx",&sna->type.txx))sna->type.txx=0;
if(!getparint("sna_type_tzz",&sna->type.tzz))sna->type.tzz=0;
if(!getparint("sna_type_txz",&sna->type.txz))sna->type.txz=0;
if(!getparint("sna_type_pp",&sna->type.pp))sna->type.pp=0;
if(!getparint("sna_type_ss",&sna->type.ss))sna->type.ss=0;
}else{
sna->type.txx=0;sna->type.tzz=0;sna->type.txz=0;
sna->type.pp=0;sna->type.ss=0;
}
/********************************/
/* Process Snapshot Information */
/********************************/
if(sna->tsnap2>=sna->tsnap1){
if(sna->vxtime<0||sna->vxtime>1) verr("snap_vxtime is not 0 or 1");
if(sna->vztime<0||sna->vztime>1) verr("snap_vztime is not 0 or 1");
if(sna->vxshift<0||sna->vxshift>1) verr("snap_vxshift is not 0 or 1");
if(sna->vzshift<0||sna->vzshift>1) verr("snap_vzshift is not 0 or 1");
if(sna->withbnd<0||sna->withbnd>1) verr("snap_snapwithbnd is not 0 or 1");
if(sna->forw<0||sna->forw>1) verr("snap_forw is not 0 or 1");
if(sna->back<0||sna->back>1) verr("snap_forw is not 0 or 1");
if(sna->beam<0||sna->beam>1) verr("snap_beam is not 0 or 1");
sna->tracl=1;
sna->isnap=0;
sna->ntr=0;
// t
sna->dtskip=(size_t)(sna->dt/mod->dt+0.5);
if(!sna->dtskip)sna->dtskip=1;
sna->dt=sna->dtskip*mod->dt;
sna->t1=(size_t)(sna->tsnap1/mod->dt+0.5);
sna->tsnap1=sna->t1*mod->dt;
sna->t2=(size_t)ceil(sna->tsnap2/mod->dt);
sna->t2-=(sna->t2-sna->t1)%sna->dtskip;
sna->nsnap=(sna->t2-sna->t1)/sna->dtskip+1;
sna->tsnap2=sna->t2*mod->dt;
// x
sna->dxskip=(size_t)(sna->dx/mod->dx);
sna->dx=sna->dxskip*mod->dx;
// z
sna->dzskip=(size_t)(sna->dz/mod->dz);
sna->dz=sna->dzskip*mod->dz;
if(sna->withbnd){
// x
if(sna->xsnap1<mod->origx-(mod->ioPx+1)*mod->dx||sna->xsnap1>mod->xmax+(mod->nax-1-mod->iePx)*mod->dx) verr("xsnap1 lies outside the modelling grid.");
if(sna->xsnap2<mod->origx-(mod->ioPx+1)*mod->dx||sna->xsnap2>mod->xmax+(mod->nax-1-mod->iePx)*mod->dx) verr("xsnap2 lies outside the modelling grid.");
sna->x1=(size_t)((sna->xsnap1-(mod->origx-mod->ioPx*mod->dx))/mod->dx+0.5);
sna->xsnap1=mod->origx+sna->x1*mod->dx-mod->ioPx*mod->dx;
sna->x2=(size_t)((sna->xsnap2-(mod->origx-mod->ioPx*mod->dx))/mod->dx+0.5);
sna->x2+=(sna->x2-sna->x1)%sna->dxskip;
if(sna->x2>mod->nax) sna->x2-=sna->dxskip;
sna->xsnap2=mod->origx+sna->x2*mod->dx-mod->ioPx*mod->dx;
sna->nx=(sna->x2-sna->x1)/sna->dxskip+1;
//z
if(sna->zsnap1<mod->origz-(mod->ioPz+1)*mod->dz||sna->zsnap1>mod->zmax+(mod->naz-1-mod->iePz)*mod->dz) verr("xsnap1 lies outside the modelling grid."); if(sna->zsnap2<mod->origz-(mod->ioPz+1)*mod->dz||sna->zsnap2>mod->zmax+(mod->naz-1-mod->iePz)*mod->dz) verr("xsnap2 lies outside the modelling grid.");
sna->z1=(size_t)((sna->zsnap1-(mod->origz-mod->ioPz*mod->dz))/mod->dz+0.5);
sna->zsnap1=mod->origz+sna->z1*mod->dz-mod->ioPz*mod->dz;
sna->z2=(size_t)((sna->zsnap2-(mod->origz-mod->ioPz*mod->dz))/mod->dz+0.5);
sna->z2+=(sna->z2-sna->z1)%sna->dzskip;
if(sna->z2>mod->naz) sna->z2-=sna->dzskip;
sna->zsnap2=mod->origz+sna->z2*mod->dz-mod->ioPx*mod->dx;
sna->nz=(sna->z2-sna->z1)/sna->dzskip+1;
}else{
// x
if(sna->xsnap1<mod->origx||sna->xsnap1>mod->xmax) verr("xsnap1 lies outside the model.");
if(sna->xsnap2<mod->origx||sna->xsnap2>mod->xmax) verr("xsnap2 lies outside the model.");
sna->x1=(size_t)((sna->xsnap1-mod->origx)/mod->dx+0.5);
sna->xsnap1=mod->origx+sna->x1*mod->dx;
sna->x2=(size_t)((sna->xsnap2-mod->origx)/mod->dx+0.5);
sna->x2+=(sna->x2-sna->x1)%sna->dxskip;
if(sna->x2>mod->nx) sna->x2-=sna->dxskip;
sna->xsnap2=mod->origx+sna->x2*mod->dx;
sna->nx=(sna->x2-sna->x1)/sna->dxskip+1;
sna->x1+=mod->ioPx;sna->x2+=mod->ioPx;
// z
if(sna->zsnap1<mod->origz||sna->zsnap1>mod->zmax) verr("zsnap1 lies outside the model.");
if(sna->zsnap2<mod->origz||sna->zsnap2>mod->zmax) verr("zsnap2 lies outside the model.");
sna->z1=(size_t)((sna->zsnap1-mod->origz)/mod->dz+0.5);
sna->zsnap1=mod->origz+sna->z1*mod->dz;
sna->z2=(size_t)((sna->zsnap2-mod->origz)/mod->dz+0.5);
sna->z2+=(sna->z2-sna->z1)%sna->dzskip;
if(sna->z2>mod->nz) sna->z2-=sna->dzskip;
sna->zsnap2=mod->origz+sna->z2*mod->dz;
sna->nz=(sna->z2-sna->z1)/sna->dzskip+1;
sna->z1+=mod->ioPz;sna->z2+=mod->ioPz;
}
/* Decomposed Snapshots? */
sna->decomp=0;
if(sna->type.pu){sna->decomp=1;decomp->decomp=1;decomp->pu=1;}
if(sna->type.pd){sna->decomp=1;decomp->decomp=1;decomp->pd=1;}
if(sna->type.pl){sna->decomp=1;decomp->decomp=1;decomp->pl=1;}
if(sna->type.pr){sna->decomp=1;decomp->decomp=1;decomp->pr=1;}
if(sna->type.pn){sna->decomp=1;decomp->decomp=2;decomp->pn=1;}
}else{
sna->nsnap=0;
sna->t1=1;
sna->t2=0;
sna->dtskip=1;
}
/****************************/
/* Process Beam Information */
/****************************/
// if(sna->beam){
// sna->skipdx=MAX(1,NINT(sna->dxsnap/mod->dx));
// sna->skipdz=MAX(1,NINT(sna->dzsnap/mod->dz));
// sna->x1=NINT((MIN(MAX(mod->origx,sna->xsnap1),mod->xmax)-mod->origx)/mod->dx);
// sna->x2=NINT((MIN(MAX(mod->origx,sna->xsnap2),mod->xmax)-mod->origx)/mod->dx);
// sna->z1=NINT((MIN(MAX(mod->origz,sna->zsnap1),mod->zmax)-mod->origz)/mod->dz);
// sna->z2=NINT((MIN(MAX(mod->origz,sna->zsnap2),mod->zmax)-mod->origz)/mod->dz);
// sna->dxsnap=mod->dx*sna->skipdx;
// sna->dzsnap=mod->dz*sna->skipdz;
// sna->nx=1+(((sna->x2-sna->x1))/sna->skipdx);
// sna->nz=1+(((sna->z2-sna->z1))/sna->skipdz);
// }
/***************************************/
/* Process Migration Image Information */
/***************************************/
if(mig->mode<0||mig->mode>5){
vwarn("Unknown migration mode (%d) expected 0-5.",mig->mode);
vwarn("Setting migration mode to conventional (1).");
mig->mode=1;
} if(mig->orient<0||(mig->orient>3&&mig->orient<9991)||mig->orient>9993){ //TODO: Test this properly and possibly remove
vwarn("Unknown migration orientation (%d) expected 0-3.",mig->orient);
vwarn("Setting migration orientation to up-down (1).");
mig->orient=1;
}
if(!mig->mode&&!(rec->rec||rec->file_loc||sna->tsnap2>=sna->tsnap1)) verr("Specified modelling mode (mig. mode=0), but no receivers or snapshots!");
if(mig->mode&&!rcv->file_src&&!rec->rec) verr("Receiver File Or Recording Grid Required For Backward Propagation!");
if(mig->mode==2&&mig->orient==3){
vwarn("Normal orientation not yet implemented for Poynting decomposition!");
vwarn("Setting orientation to 0.");
mig->orient=0;
}
if(mig->mode==4&&mig->orient!=1){
vwarn("Only up-down imaging currently implemented for plane-wave decompostion");
vwarn("Setting orientation to up-down (1) imaging");
mig->orient=1;
}
if(mig->mode==5&&(mig->orient<1||mig->orient>2)){
vwarn("For Hilbert transform based imaging only up-down (1) or left-right (2)");
vwarn("imaging possible.");
vwarn("Setting orientation to up-down (1) imaging");
mig->orient=1;
}
mig->skipdt=(size_t)(mig->dt/mod->dt+0.5);
if(!mig->skipdt){mig->skipdt=1;vwarn("Migration time step (%f) less than modelling time step (%f) increasing to modelling time step!",mig->dt,mod->dt);}
mig->dt=mig->skipdt*mod->dt;
mig->skipdx=(size_t)(mig->dx/mod->dx+0.5);
mig->dx=mig->skipdx*mod->dx;
// mig->nx=mig->nx/mig->skipdx+1;
mig->skipdz=(size_t)(mig->dz/mod->dz+0.5);
mig->dz=mig->skipdz*mod->dz;
// mig->nz=mig->nz/mig->skipdz+1;
if(mig->xmig1<mod->origx||mig->xmig1>mod->xmax) vwarn("xmig1 (%f) lies outside the model.",mig->xmig1);
if(mig->xmig2<mod->origx||mig->xmig2>mod->xmax) vwarn("xmig2 (%f) lies outside the model.",mig->xmig2);
if(mig->zmig1<mod->origz||mig->zmig1>mod->xmax) vwarn("zmig1 (%f) lies outside the model.",mig->zmig1);
if(mig->zmig2<mod->origz||mig->zmig2>mod->xmax) vwarn("zmig2 (%f) lies outside the model.",mig->zmig2);
mig->x1=(size_t)((MIN(MAX(mod->origx,mig->xmig1),mod->xmax)-mod->origx)/mod->dx+0.5);
mig->x2=(size_t)((MIN(MAX(mod->origx,mig->xmig2),mod->xmax)-mod->origx)/mod->dx+0.5);
mig->z1=(size_t)((MIN(MAX(mod->origz,mig->zmig1),mod->zmax)-mod->origz)/mod->dz+0.5);
mig->z2=(size_t)((MIN(MAX(mod->origz,mig->zmig2),mod->zmax)-mod->origz)/mod->dz+0.5);
if(mig->skipdt>1)mig->nt=mod->nt/mig->skipdt+1;else mig->nt=mod->nt;
mig->nx=(mig->x2-mig->x1)/mig->skipdx+1;
mig->nz=(mig->z2-mig->z1)/mig->skipdz+1;
mig->sizem=mig->nx*mig->nz;
if(mig->nx==1) vwarn("Horizontal migration image size is only 1.");
if(mig->nz==1) vwarn("Vertical migration image size is only 1.");
mig->ntr=0;
if(mig->mode==3){ //Decomposition Mode
switch(mig->orient){
case 1: //Up-Down
decomp->decomp=1; //Turn Snapshot-Decomposition On
decomp->pu=1; //Up-Going Decomposition
decomp->pd=1; //Down-Going Decomposition
break;
case 2: //Left-Right
decomp->decomp=1; //Turn Snapshot-Decomposition On
decomp->pl=1; //Left-Going Decomposition
decomp->pr=1; //Right-Going Decomposition
break;
case 3: //Normal
decomp->decomp=2; //Turn Snapshot-Decomposition On - 2 for directional decomposition
decomp->pn=1; //Normal-Going Decomposition
break;
case 9991: // Up-Down -- Individual Components
decomp->decomp=1;
decomp->pu=1; //Up-Going Decomposition
decomp->pd=1; //Down-Going Decomposition
break;
case 9992: //Left-Right -- Individual Components
//TODO: Implement break;
case 9993: //Normal -- Individual Components
//TODO: Implement
break;
}
}
if(!mig->tap) mig->tap=MIN(MIN(mig->nx,mig->nt),mig->nz)/10;
/**************************************/
/* Initialize Random Number Generator */
/**************************************/
srand((unsigned)time(NULL)); //Random Number Seed
/***************************************/
/* Prepare Finite Difference Operators */
/***************************************/
prepareFDOperators(mod,bnd,decomp);
/******************/
/* Initialize PML */
/******************/
if(bnd->pml) initPML(mod,bnd,*verbose);
/**********************************/
/* Dispersion & Stability Factors */
/**********************************/
if(mod->iorder==2){dispfactor=10;stabfactor=1.0/sqrt(2.0);}else{dispfactor=5;stabfactor=0.606;}
/*************************************/
/* Process Decomposition Information */
/*************************************/
if(decomp->wavFilt){
decomp->kl*=mod->dx*dispfactor/12; //Determine start of filter
decomp->kh*=mod->dx*dispfactor/12; //Determine end of filter
}
/***************************/
/* Initialize Compression? */
/***************************/
if(!getparint("compress",&(mig->compress))) mig->compress=0; //Default: No Compression
else if(mig->compress){
initializeCompressionField(mig); //Initialize compresion information
}
if(*verbose){
vmess("*******************************************");
vmess("****** Finite Difference Elastic RTM ******");
if(*verbose>3) vmess("******* RTM: Reverse Time Migration *******");
vmess("*********** General Information ***********");
vmess("*******************************************");
if(mod->ischeme==1) vmess("Acoustic staggered grid, pressure/velocity");
if(mod->ischeme==2) vmess("Visco-Acoustic staggered grid, pressure/velocity");
if(mod->ischeme==3) vmess("Elastic staggered grid, stress/velocity");
if(mod->ischeme==4) vmess("Visco-Elastic staggered grid, stress/velocity");
vmess("dt= %dus dt= %fs(%e)",mod->dtus,mod->dt,mod->dt);
// vmess("tmod = %f",mod->tmod);
// vmess("ntsam= %d dt= %f(%e)",mod->nt, mod->dt, mod->dt);
vmess("*******************************************");
vmess("************ Model Information ************");
vmess("*******************************************");
vmess("P-wave velocity file: %s",mod->file_cp);
if(mod->ischeme>2) vmess("S-wave velocity file: %s",mod->file_cs);
if(mod->file_den) vmess("Density file: %s",mod->file_den);
if(mod->file_imp)vmess("* Impedance file: %s",mod->file_imp);
vmess("Model Dimensions:");
vmess("nz = %8d nx = %8d",mod->nz,mod->nx);
vmess("dz = %8.4f dx = %8.4f",mod->dz,mod->dx);
vmess("zmin = %8.4f zmax = %8.4f",mod->origz,mod->zmax);
vmess("xmin = %8.4f xmax = %8.4f",mod->origx,mod->xmax);
vmess("min(cp) = %9.3f max(cp) = %9.3f",mod->cp_min,mod->cp_max); if(mod->ischeme>2)vmess("min(cs) = %9.3f max(cs) = %9.3f",mod->cs_min,mod->cs_max);
vmess("min(ro) = %9.3f max(ro) = %9.3f",mod->rho_min,mod->rho_max);
vmess("*******************************************");
vmess("********* Dispersion & Stability **********");
vmess("*******************************************");
vmess("Dispersion criterion is %3d points per wavelength: ", NINT(dispfactor));
vmess(" ====> wavelength > %f m [dx*disp]",mod->dx*dispfactor);
vmess("The maximum frequency in source wavelet must be:");
vmess(" ====> frequency < %f Hz. [Cmin/dx*disp]",cmin/(mod->dx*dispfactor));
vmess("Stability criterion for current settings: ");
vmess(" ====> Cp < %f m/s [dx*disp/dt]",mod->dx*stabfactor/mod->dt);
vmess(" Optimal discretisation for current model:");
vmess(" With maximum velocity = %f dt <= %e", cmax,mod->dx*stabfactor/cmax);
vmess(" With grid spacing = %f fmax <= %e", mod->dx, cmin/(mod->dx*dispfactor));
vmess("*******************************************");
vmess("************** Boundary Info **************");
vmess("*******************************************");
vmess("*** 1=free 2=pml 3=rigid 4=tapered ***");
vmess("Top boundary : %d",bnd->top);
vmess("Left boundary : %d",bnd->lef);
vmess("Right boundary : %d",bnd->rig);
vmess("Bottom boundary : %d",bnd->bot);
if(bnd->npml!=0) vmess("PML length = %d points",bnd->npml);
if(bnd->ntapo!=0) vmess("Taper length = %d points",bnd->ntapo);
if(bnd->npml!=0&&bnd->ntapo!=0&&bnd->npml!=bnd->ntapo){
vmess("Both PML & tapered boundaries are used.");
vmess("Hence the larger boundary will be used.");
vmess("Boundary size nbnd=%d");
}
if(*verbose>2){
vmess("*******************************************");
vmess("************** Modeling Grid **************");
vmess("*******************************************");
if(mod->ischeme==1) vmess("Acoustic Staggered Grid, Pressure/Velocity");
else if(mod->ischeme==2) vmess("Visco-Acoustic Staggered Grid, Pressure/Velocity");
else if(mod->ischeme==3) vmess("Elastic Staggered Grid, Stress/Velocity");
else if(mod->ischeme==4) vmess("Visco-Elastic Staggered Grid, Stress/Velocity");
vmess("Corner Grid-Point Indices:");
printf(" %s: %c┌─────────────────────────────────────────┐%c\n",xargv[0],14,15);
printf(" %s: %c│E1 EDGE │%c\n",xargv[0],14,15);
printf(" %s: %c│ ┌─────────────────────────────────┐ │%c\n",xargv[0],14,15);
printf(" %s: %c│ │B1 BOUNDARY │ │%c\n",xargv[0],14,15);
printf(" %s: %c│ │ ┌─────────────────────────┐ │ │%c\n",xargv[0],14,15);
printf(" %s: %c│ E │ │M1 │ │ E │%c\n",xargv[0],14,15);
printf(" %s: %c│ D │ B │ │ B │ D │%c\n",xargv[0],14,15);
printf(" %s: %c│ G │ N │ Modeling Domain │ N │ G │%c\n",xargv[0],14,15);
printf(" %s: %c│ E │ D │ │ D │ E │%c\n",xargv[0],14,15);
printf(" %s: %c│ │ │ M2│ │ │%c\n",xargv[0],14,15);
printf(" %s: %c│ │ └─────────────────────────┘ │ │%c\n",xargv[0],14,15);
printf(" %s: %c│ │ BOUNDARY B2│ │%c\n",xargv[0],14,15);
printf(" %s: %c│ └─────────────────────────────────┘ │%c\n",xargv[0],14,15);
printf(" %s: %c│ EDGE E2│%c\n",xargv[0],14,15);
printf(" %s: %c└─────────────────────────────────────────┘%c\n",xargv[0],14,15);
if(mod->ischeme<3) vmess("Pressure Grid:");else vmess("Txx/Tzz Grid:");
vmess("E1:(0,0) E2:(%d,%d)",mod->nax-2,mod->naz-2);
vmess("B1:(%d,%d) B2:(%d,%d)",mod->ioPxb,mod->ioPzb,mod->iePxb-1,mod->iePzb-1);
vmess("M1:(%d,%d) M2:(%d,%d)",mod->ioPx ,mod->ioPz, mod->iePx-1 ,mod->iePz-1 );
if(mod->ischeme>2){
vmess("Txz Grid:");
vmess("E1:(0,0) E2:(%d,%d)",mod->nax-1,mod->naz-1);
vmess("M1:(%d,%d) M2:(%d,%d)",mod->ioTxb,mod->ioTzb,mod->ieTxb-1,mod->ieTzb-1);
vmess("M1:(%d,%d) M2:(%d,%d)",mod->ioTx ,mod->ioTz, mod->ieTx-1 ,mod->ieTz-1 );
}
vmess("Vx Grid:");
vmess("E1:(0,0) E2:(%d,%d)",mod->nax-1,mod->naz-2);
vmess("B1:(%d,%d) B2:(%d,%d)",mod->ioXxb,mod->ioXzb,mod->ieXxb-1,mod->ieXzb-1);
vmess("M1:(%d,%d) m2:(%d,%d)",mod->ioXx ,mod->ioXz, mod->ieXx-1 ,mod->ieXz-1 );
vmess("Vz Grid:");
vmess("E1:(0,0) E2:(%d,%d)",mod->nax-2,mod->naz-1);
vmess("B1:(%d,%d) B2:(%d,%d)",mod->ioZxb,mod->ioZzb,mod->ieZxb-1,mod->ieZzb-1); vmess("M1:(%d,%d) M2:(%d,%d)",mod->ioZx ,mod->ioZz, mod->ieZx-1 ,mod->ieZz-1 );
}
vmess("*******************************************");
vmess("********** Source Wavefield Info **********");
vmess("*******************************************");
vmess("Source wavefield file: %s",src->file_src);
// vmess("Number of sources defined: %d",src->nsrc);
vmess("Sources where applicable were moved to the");
vmess("nearest grid point.");
// if(src->nsrc<=10||*verbose>3){
// vmess("Source Type:");
// vmess("1=P 2=Txz 3=Tzz 4=Txx 5=S-pot 6=Fx 7=Fz 8=P-pot");
// vmess("Source Orientation:");
// vmess("1=Monopole 2=Vertical Dipole -/+ 3=Horizontal Dipole -/+");
// for(ix=0;ix<src->nsrc;ix++) vmess("Source %d at (%f,%f) of type %d and orientation %d.",ix+1,src->x[ix],src->z[ix],src->typ[ix],src->orient[ix]);
// }
vmess("*******************************************");
vmess("********* Receiver Wavefield Info *********");
vmess("*******************************************");
if(!mig->mode&&rec->rec&&rec->write){ //Modelling Data?
vmess("Receiver wavefield file: %s",rcv->file_src);
if(rec->top||rec->left||rec->right||rec->bottom){
vmess("Receivers on the following boundaries:");
fprintf(stderr," %s:",xargv[0]);
if(rec->top) fprintf(stderr," top");
if(rec->left) fprintf(stderr," left");
if(rec->right) fprintf(stderr," right");
if(rec->bottom)fprintf(stderr," bottom");
fprintf(stderr,"\n");
if(rec->file_loc) vmess("With additional receivers defined in: %s\n",rec->file_loc);
}else vmess("Receiver locations defined in: %s",rec->file_loc);
}else{
if(rcv->file_src) vmess("Receiver wavefield file: %s",rcv->file_src);
else{
vmess("Receivers on the following boundaries:");
fprintf(stderr," %s:",xargv[0]);
if(rec->top) fprintf(stderr," top");
if(rec->left) fprintf(stderr," left");
if(rec->right) fprintf(stderr," right");
if(rec->bottom)fprintf(stderr," bottom");
fprintf(stderr,"\n");
if(rec->file_loc) vmess("With additional receivers defined in: %s\n",rec->file_loc);
}
}
if(rec->rec) vmess("Number of receivers defined: %d",rcv->nsrc);
vmess("Receivers where applicable were moved to");
vmess("the nearest grid point.");
if(rec->rec&&(rcv->nsrc<=10||*verbose>3)){
vmess("Receiver Type:");
vmess("1=P 2=Txz 3=Tzz 4=Txx 5=S-pot 6=Fx 7=Fz 8=P-pot");
vmess("Receiver Orientation:");
vmess("1=Monopole 2=Vertical Dipole -/+ 3=Horizontal Dipole -/+");
for(ix=0;ix<rcv->nsrc;ix++) vmess("Receiver %d at (%f,%f) [%d,%d] of type %d and orientation %d.",ix+1,rcv->x[ix],rcv->z[ix],rcv->xi[ix],rcv->zi[ix],rcv->typ[ix],rcv->orient[ix]);
}
if(rec->write) vmess("Receiver data will be written to disk.");
if(sna->nsnap>0){
vmess("*******************************************");
vmess("************* snap shot info **************");
vmess("*******************************************");
vmess("Snapshot Base File: %s",sna->file_snap);
vmess("tsnap1=%f tsnap2=%f",sna->tsnap1,sna->tsnap2);
vmess("dtsnap=%f Nsnap =%d",sna->dt,sna->nsnap);
vmess("nxsnap=%d nzsnap=%d",sna->nx,sna->nz);
vmess("dxsnap=%f dzsnap=%f",sna->dx,sna->dz);
vmess("xmin =%f xmax =%f",sna->xsnap1,sna->xsnap2);
vmess("zmin =%f zmax =%f",sna->zsnap1,sna->zsnap2);
if(sna->vxtime) vmess("vx snapshot time : t+0.5*dt");
else vmess("vx snapshot time : t-0.5*dt ");
if(sna->vztime) vmess("vz snapshot time : t+0.5*dt");
else vmess("vz snapshot time : t-0.5*dt "); fprintf(stderr," %s: Snapshot types : ",xargv[0]);
if(sna->type.vz) fprintf(stderr,"Vz ");
if(sna->type.vx) fprintf(stderr,"Vx ");
if(mod->ischeme<3){ //Acoustic
if(sna->type.p) fprintf(stderr,"P ");
if(sna->type.pu) fprintf(stderr,"Pu ");
if(sna->type.pd) fprintf(stderr,"Pd ");
if(sna->type.pl) fprintf(stderr,"Pl ");
if(sna->type.pr) fprintf(stderr,"Pr ");
if(sna->type.pn) fprintf(stderr,"Pn ");
}else{ //Elastic
if(sna->type.txx) fprintf(stderr,"Txx ");
if(sna->type.tzz) fprintf(stderr,"Tzz ");
if(sna->type.txz) fprintf(stderr,"Txz ");
if(sna->type.pp) fprintf(stderr,"P ");
if(sna->type.ss) fprintf(stderr,"S ");
}
fprintf(stderr,"\n");
}else{
vmess("*******************************************");
vmess("*************** No Snapshots **************");
vmess("*******************************************");
}
if(sna->beam){
vmess("*******************************************");
vmess("**************** beam info ****************");
vmess("*******************************************");
vmess("nzsnap=%d nxsnap=%d",sna->nz,sna->nx);
vmess("dzsnap=%f dxsnap=%f",sna->dz,sna->dx);
vmess("zmin =%f zmax =%f",sna->zsnap1,sna->zsnap2);
vmess("xmin =%f xmax =%f",sna->xsnap1,sna->xsnap2);
fprintf(stderr," %s: Beam types : ",xargv[0]);
if(sna->type.vz) fprintf(stderr,"Vz ");
if(sna->type.vx) fprintf(stderr,"Vx ");
if(mod->ischeme<3){
if(sna->type.p ) fprintf(stderr,"P ");
if(sna->type.pu) fprintf(stderr,"Pu");
if(sna->type.pd) fprintf(stderr,"Pd");
if(sna->type.pl) fprintf(stderr,"Pl");
if(sna->type.pr) fprintf(stderr,"Pr");
if(sna->type.pn) fprintf(stderr,"Pn");
}else if(mod->ischeme<5){
if(sna->type.txx) fprintf(stderr,"Txx ");
if(sna->type.tzz) fprintf(stderr,"Tzz ");
if(sna->type.txz) fprintf(stderr,"Txz ");
if(sna->type.pp) fprintf(stderr,"P ");
if(sna->type.ss) fprintf(stderr,"S ");
}
fprintf(stderr,"\n");
}else{
vmess("*******************************************");
vmess("**************** No Beams *****************");
vmess("*******************************************");
}
vmess("*******************************************");
vmess("******* Migration Image Information *******");
vmess("*******************************************");
if(mig->mode){
vmess("Migration image file: %s",mig->file_mig);
switch(mig->mode){
case 1:
vmess("Conventional RTM imaging condition.");
break;
case 2:
vmess("Poynting RTM imaging condition.");
break;
case 3:
vmess("Decomposition RTM imaging condition.");
break;
case 4: vmess("Plane-Wave RTM imaging condition.");
break;
case 5:
vmess("Hilbert Transform RTM imaging condition.");
break;
}
switch(mig->orient){
case 0:
vmess("Imaging wavefields not travelling in the same direction.");
break;
case 1:
vmess("Up-Down Imaging.");
break;
case 2:
vmess("Left-Right Imaging.");
break;
case 3:
vmess("Surface Normal Imaging.");
break;
}
if(mig->backscatter) vmess("Imaging Backscattered Source Receiver Wavefields.");
else vmess("Imaging Opposing Source Receiver Wavefields.");
vmess("Migration Image Information:");
vmess("nzmig=%d nxmig=%d",mig->nz,mig->nx);
vmess("dzmig=%f dxmig=%f",mig->dz,mig->dx);
vmess("xmin =%f xmax =%f",mig->xmig1,mig->xmig2);
vmess("zmin =%f zmax =%f",mig->zmig1,mig->zmig2);
vmess("Migration Snapshot Information:");
vmess("dtmig=%f",mig->dt);
// vmess("ntmig=%d dtmig=%f",mig->nt,mig->dt);
if(mig->compress) vmess("Compressing Migration Snapshots!");
}else{
vmess("No Migration. Only Modelling.");
}
if(decomp->decomp){
vmess("*******************************************");
vmess("*** Directional Wavefield Decomposition ***");
vmess("*******************************************");
vmess("* Tikhonov Regularization: %9.8e *",decomp->reg);
if(decomp->med==3){ // Check median filter type
vmess("* Median Filter Order: 3x3 *");
}else if(decomp->med==5){
vmess("* Median Filter Order: 5x5 *");
}else if(decomp->med==1){
vmess("* No median filter will be applied. *");
decomp->med=0;
}else if(decomp->med!=0){
vmess("* Median Filter Order: UNKOWN *");
vmess("* No median filter will be applied. *");
decomp->med=0;
}else{
vmess("* No median filter will be applied. *");
}
if(decomp->wavFilt){
vmess("*******************************************");
vmess("* Operator will be wavenumber filtered. *");
vmess("* Cosine Squared Filter *");
vmess("* from %09.4f to %09.4f. *",decomp->kl,decomp->kh);
if(*verbose>2) printWisdomPath();
}
if(decomp->decomp==2){
vmess("*******************************************");
vmess("* Directional Decomposition Information *");
vmess("* Preferential Direction Information *");
vmess("* px=%11.9f, pz=%11.9f *",decomp->px,decomp->pz);
if(mod->file_dd){
vmess("* Decomposition Direction File Name: %s",mod->file_dd);
if(!decomp->writeDD) vmess("* Loading Decomp. Direction From File. *");
else{ vmess("* Estimating Decomp. Dir. From Impedance. *");
vmess("* Writing Decomposition Direction To File.*"); if(mod->mavgi==3) vmess("* Impedance Average Filter Order: 3x3 *");
else if(mod->mavgi==5) vmess("* Impedance Average Filter Order: 5x5 *");
else if(mod->mavgi==7) vmess("* Impedance Average Filter Order: 7x7 *");
else if(mod->mavgi==9) vmess("* Impedance Average Filter Order: 9x9 *");
else vmess("* No Impedance Average Filter *");
if(decomp->mavgn==3) vmess("* Preferential Average Filter Order: 3x3 *");
else if(decomp->mavgn==5)vmess("* Preferential Average Filter Order: 5x5 *");
else if(decomp->mavgn==7)vmess("* Preferential Average Filter Order: 7x7 *");
else if(decomp->mavgn==9)vmess("* Preferential Average Filter Order: 9x9 *");
else vmess("* No Preferential Average Filter *");
if(decomp->mavgo==3) vmess("* Orthogonal Average Filter Order: 3x3 *");
else if(decomp->mavgo==5)vmess("* Orthogonal Average Filter Order: 5x5 *");
else if(decomp->mavgo==7)vmess("* Orthogonal Average Filter Order: 7x7 *");
else if(decomp->mavgo==9)vmess("* Orthogonal Average Filter Order: 9x9 *");
else vmess("* No Orthogonal Average Filter *");
if(decomp->mavga==3) vmess("* Dir. Angle Average Filter Order: 3x3 *");
else if(decomp->mavga==5)vmess("* Dir. Angle Average Filter Order: 5x5 *");
else if(decomp->mavga==7)vmess("* Dir. Angle Average Filter Order: 7x7 *");
else if(decomp->mavga==9)vmess("* Dir. Angle Average Filter Order: 9x9 *");
else vmess("* No Dir. Angle Average Filter *");
}
}else{
vmess("* Estimating Decomp. Dir. From Impedance. *");
if(mod->mavgi==3) vmess("* Impedance Average Filter Order: 3x3 *");
else if(mod->mavgi==5) vmess("* Impedance Average Filter Order: 5x5 *");
else if(mod->mavgi==7) vmess("* Impedance Average Filter Order: 7x7 *");
else if(mod->mavgi==9) vmess("* Impedance Average Filter Order: 9x9 *");
else vmess("* No Impedance Average Filter *");
if(decomp->mavgn==3) vmess("* Preferential Average Filter Order: 3x3 *");
else if(decomp->mavgn==5)vmess("* Preferential Average Filter Order: 5x5 *");
else if(decomp->mavgn==7)vmess("* Preferential Average Filter Order: 7x7 *");
else if(decomp->mavgn==9)vmess("* Preferential Average Filter Order: 9x9 *");
else vmess("* No Preferential Average Filter *");
if(decomp->mavgo==3) vmess("* Orthogonal Average Filter Order: 3x3 *");
else if(decomp->mavgo==5)vmess("* Orthogonal Average Filter Order: 5x5 *");
else if(decomp->mavgo==7)vmess("* Orthogonal Average Filter Order: 7x7 *");
else if(decomp->mavgo==9)vmess("* Orthogonal Average Filter Order: 9x9 *");
else vmess("* No Orthogonal Average Filter *");
if(decomp->mavga==3) vmess("* Dir. Angle Average Filter Order: 3x3 *");
else if(decomp->mavga==5)vmess("* Dir. Angle Average Filter Order: 5x5 *");
else if(decomp->mavga==7)vmess("* Dir. Angle Average Filter Order: 7x7 *");
else if(decomp->mavga==9)vmess("* Dir. Angle Average Filter Order: 9x9 *");
else vmess("* No Dir. Angle Average Filter *");
}
}
vmess("*******************************************");
// vmess("* The following decomposed data will be *");
// vmess("* be written out: *");
// // Snapshots
// if(decomp->sPu)vmess("* Up-going pressure snapshots *");
// if(decomp->sPd)vmess("* Down-going pressure snapshots *");
// if(decomp->sPl)vmess("* Left-going pressure snapshots *");
// if(decomp->sPr)vmess("* Right-going pressure snapshots *");
// if(decomp->sPn)vmess("* Normal-going pressure snapshots *");
// if(decomp->sVxu)vmess("* Up-going vx snapshots *");
// if(decomp->sVxd)vmess("* Down-going vx snapshots *");
// if(decomp->sVxl)vmess("* Left-going vx snapshots *");
// if(decomp->sVxr)vmess("* Right-going vx snapshots *");
// if(decomp->sVxn)vmess("* Normal-going vx snapshots *");
// if(decomp->sVzu)vmess("* Up-going vz snapshots *");
// if(decomp->sVzd)vmess("* Down-going vz snapshots *");
// if(decomp->sVzl)vmess("* Left-going vz snapshots *");
// if(decomp->sVzr)vmess("* Right-going vz snapshots *");
// if(decomp->sVzn)vmess("* Normal-going vz snapshots *");
// // Receiver Data
// if(decomp->rPu)vmess("* Up-going pressure receivers *");
// if(decomp->rPd)vmess("* Down-going pressure receivers *");
// if(decomp->rPl)vmess("* Left-going pressure receivers *");
// if(decomp->rPr)vmess("* Right-going pressure receivers *");
// if(decomp->rPn)vmess("* Normal-going pressure receivers *");// if(decomp->rVxu)vmess("* Up-going vx receivers *");
// if(decomp->rVxd)vmess("* Down-going vx receivers *");
// if(decomp->rVxl)vmess("* Left-going vx receivers *");
// if(decomp->rVxr)vmess("* Right-going vx receivers *");
// if(decomp->rVxn)vmess("* Normal-going vx receivers *");
// if(decomp->rVzu)vmess("* Up-going vz receivers *");
// if(decomp->rVzd)vmess("* Down-going vz receivers *");
// if(decomp->rVzl)vmess("* Left-going vz receivers *");
// if(decomp->rVzr)vmess("* Right-going vz receivers *");
// if(decomp->rVzn)vmess("* Normal-going vz receivers *");
}else{
vmess("*******************************************");
vmess("** No Directional Wavefield Decomposition *");
vmess("*******************************************");
}
}
/****************************/
/* Correct For Free Surface */
/****************************/
if(bnd->top==1)mod->ioPz++;
if(bnd->lef==1)mod->ioPx++;
if(bnd->rig==1)mod->iePx--;
if(bnd->bot==1)mod->iePz--;
/****************************/
/* Correct For Rigid Surface */
/****************************/
if(bnd->top==3)mod->ioXz++;
if(bnd->lef==3)mod->ioZx++;
if(bnd->rig==3)mod->ieZx--;
if(bnd->bot==3)mod->ieXz--;
return(0);
}