...
 
Commits (11)
......@@ -4,4 +4,5 @@
*.backup
*~*
*.orig
*.o
manual_IFOS3D.pdf
......@@ -24,8 +24,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
......@@ -55,7 +54,7 @@ float *** taus, float *** taup, float * eta){
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
......@@ -86,12 +85,12 @@ float *** taus, float *** taup, float * eta){
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
if(L){
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......@@ -118,12 +117,12 @@ float *** taus, float *** taup, float * eta){
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
if(L){
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......
......@@ -5,8 +5,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
......@@ -36,7 +35,7 @@ float *** taus, float *** taup, float * eta){
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
......@@ -70,12 +69,13 @@ float *** taus, float *** taup, float * eta){
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
if(L){
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......@@ -92,7 +92,7 @@ float *** taus, float *** taup, float * eta){
sprintf(modfile,"%s.mod",MFILE);
writemod(modfile,rho,3);
writemod(modfile,mod,3);
MPI_Barrier(MPI_COMM_WORLD);
......
......@@ -5,8 +5,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float *eta) {
void model(st_model *mod) {
/*-----------------------------------------------------------------------*/
/* extern variables */
......@@ -37,7 +36,7 @@ void model(float *** rho, float *** pi, float *** u,
for (l=1; l<=L; l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
......@@ -73,13 +72,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......
......@@ -5,8 +5,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float *eta) {
void model(st_model *mod) {
/*-----------------------------------------------------------------------*/
/* extern variables */
......@@ -38,7 +37,7 @@ void model(float *** rho, float *** pi, float *** u,
for (l=1; l<=L; l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
......@@ -74,13 +73,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......@@ -111,13 +110,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......@@ -148,13 +147,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......@@ -185,13 +184,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......@@ -222,13 +221,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......
/*
* Model defined by flnode file.
*/
#include "fd.h"
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY,NZ, NXG, NYG, NZG, POS[4], MYID;
extern float DZ;
extern FILE *FP;
/* local variables */
float muv, piv, vp, vs, rhov;
int i, j,k, ii, jj,kk, l;
FILE *flfile;
int nodes;
char cline[256];
float *fldepth, *flrho, *flvp, *flvs;
nodes=7;
fldepth=vector(1,nodes);
flrho=vector(1,nodes);
flvp=vector(1,nodes);
flvs=vector(1,nodes);
/*read FL nodes from File*/
flfile=fopen("flnodes/flnodes.rheinstetten","r");
if (flfile==NULL) err(" FL-file could not be opened !");
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
}
else l=l-1;
}
if(MYID==0){
fprintf(FP," ------------------------------------------------------------------ \n\n");
fprintf(FP," Information of FL nodes: \n\n");
fprintf(FP," \t depth \t rho \t vp \t vs \n\n");
for (l=1;l<=nodes;l++){
fprintf(FP," \t %f \t %f \t %f \t %f\n\n",fldepth[l],flrho[l],flvp[l],flvs[l]);
}
fprintf(FP," ------------------------------------------------------------------ \n\n");
}
/*-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------*/
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (k=1; k<=NZG; k++) {
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (k==1) && (MYID==0)){
fprintf(FP,"depth: %f m: double node\n",fldepth[l]);}}
else{
for (j=(int)(fldepth[l]/DZ)+1;j<=(int)(fldepth[l+1]/DZ);j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=(DZ*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
vs=(DZ*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
vs=vs*1000.0;
rhov=(DZ*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
rhov=rhov*1000.0;
muv=vs*vs*rhov;
piv=vp*vp*rhov;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=rhov;
mod->pi[jj][ii][kk]=piv;
}
}
}
}
for (j=(int)(fldepth[nodes]/DZ)+1;j<=NYG;j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=flvp[nodes]*1000.0; vs=flvs[nodes]*1000.0; rhov=flrho[nodes]*1000.0;
muv=vs*vs*rhov;
piv=vp*vp*rhov;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=rhov;
mod->pi[jj][ii][kk]=piv;
}
}
}
}
free_vector(fldepth,1,nodes);
free_vector(flrho,1,nodes);
free_vector(flvp,1,nodes);
free_vector(flvs,1,nodes);
}
......@@ -205,7 +205,7 @@ caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','normal');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at x=' num2str(sly) 'm'],'FontSize',fontsize)
title(['true ' parameter '-model. Slize at y=' num2str(sly) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
......@@ -236,7 +236,7 @@ caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','normal');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at x=' num2str(sly) 'm'],'FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at y=' num2str(sly) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
......@@ -267,7 +267,7 @@ set(gca,'ydir','normal');
%set(gca,'xdir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at x=' num2str(slz) 'm'],'FontSize',fontsize)
title(['true ' parameter '-model. Slize at z=' num2str(slz) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
......@@ -298,7 +298,7 @@ set(gca,'ydir','normal');
%set(gca,'xdir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at x=' num2str(slz) 'm'],'FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at z=' num2str(slz) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
......
# run forward simulation to obtain observed data
make clean
make install MODEL=hh_toy_true.c
make install MODEL=../genmod/hh_toy_true.c
mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/ifos3d_toy_FW.json | tee ./in_and_out/ifos3D.out
cp model/toy.vs_it0 model/toy.vs.true
......@@ -10,5 +10,5 @@ cp model/toy.rho_it0 model/toy.rho.true
# invert observed data using homogeneous starting model
make clean
make install MODEL=hh_toy_start.c
mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/ifos3d_toy.json | tee ./in_and_out/ifos3D_INV.out
make install MODEL=../genmod/hh_toy_start.c
mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/ifos3d_toy.json | tee ./in_and_out/ifos3D_INV.out
This diff is collapsed.
......@@ -8,7 +8,9 @@
# source code for model generation (elastic)
MODEL_SCR = hh_toy_true.c
MODEL_SCR = ../genmod/hh_toy_true.c
#MODEL_SCR = ../genmod/rheinstetten_elastic.c
# Compiler (LAM: CC=hcc, CRAY T3E: CC=cc, SHARCNet: mpicc)
......@@ -17,7 +19,7 @@ CC=mpicc
LFLAGS=-lm -lmpi -lcseife
CFLAGS=-Wall -O4
SFLAGS=-L./../libcseife
IFLAGS=-I./../libcseife
IFLAGS=-I./../libcseife -I.
#Hermit
#CC=cc
......@@ -44,8 +46,12 @@ IFLAGS=-I./../libcseife
# after this line, no further editing should be necessary
# --------------------------------------------------------
.c.o:
$(CC) -c $(CFLAGS) $< $(IFLAGS)
# .c.o:
# $(CC) -c $(CFLAGS) $< $(IFLAGS)
#
%.o: %.c
$(CC) $(CFLAGS) -o $@ -c $< $(IFLAGS)
# $(CC) $(CFLAGS) -MM $*.c $(IFLAGS) > $*.d
SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o)
......@@ -162,12 +168,15 @@ snapmerge: $(SNAPMERGE_OBJ)
# part_model: $(PARTMODEL_OBJ)
# $(CC) $(LFLAGS) $(PARTMODEL_OBJ) -o ../bin/partmodel
ifos3d: $(IFOS_OBJ)
$(CC) $(SFLAGS) $(IFOS_OBJ) $(LFLAGS) -o ../bin/ifos3d
ifos3d: $(IFOS_OBJ) fd.h
$(CC) $(SFLAGS) $(IFOS_OBJ) $(LFLAGS) -o ../bin/ifos3d
clean:
find . -name "*.o" -exec rm {} \;
find . -name "*.bck" -exec rm {} \;
find . -name "core" -exec rm {} \;
find ../genmod -name "*.o" -exec rm {} \;
find ../genmod -name "*.d" -exec rm {} \;
......@@ -25,10 +25,7 @@
#include "fd.h"
void av_mat(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup,
float *** uipjp, float *** ujpkp, float *** uipkp, float *** tausipjp,
float *** tausjpkp, float *** tausipkp, float *** rjp, float *** rkp, float *** rip ){
void av_mat(st_model *mod,st_model_av *mod_av){
extern int NX, NY, NZ, MYID,L,VERBOSE;
......@@ -48,20 +45,20 @@ float *** tausjpkp, float *** tausipkp, float *** rjp, float *** rkp, float
/* harmonic averaging of shear modulus */
uipjp[j][i][k]=4.0/((1.0/u[j][i][k])+(1.0/u[j][i+1][k])+(1.0/u[j+1][i+1][k])+(1.0/u[j+1][i][k]));
ujpkp[j][i][k]=4.0/((1.0/u[j][i][k])+(1.0/u[j][i][k+1])+(1.0/u[j+1][i][k+1])+(1.0/u[j+1][i][k]));
uipkp[j][i][k]=4.0/((1.0/u[j][i][k])+(1.0/u[j][i][k+1])+(1.0/u[j][i+1][k+1])+(1.0/u[j][i+1][k]));
mod_av->uipjp[j][i][k]=4.0/((1.0/mod->u[j][i][k])+(1.0/mod->u[j][i+1][k])+(1.0/mod->u[j+1][i+1][k])+(1.0/mod->u[j+1][i][k]));
mod_av->ujpkp[j][i][k]=4.0/((1.0/mod->u[j][i][k])+(1.0/mod->u[j][i][k+1])+(1.0/mod->u[j+1][i][k+1])+(1.0/mod->u[j+1][i][k]));
mod_av->uipkp[j][i][k]=4.0/((1.0/mod->u[j][i][k])+(1.0/mod->u[j][i][k+1])+(1.0/mod->u[j][i+1][k+1])+(1.0/mod->u[j][i+1][k]));
/* arithmetic averaging of TAU for S-waves and density */
if (L){
tausipjp[j][i][k]=0.25*(taus[j][i][k]+taus[j][i+1][k]+taus[j+1][i+1][k]+taus[j+1][i][k]);
tausjpkp[j][i][k]=0.25*(taus[j][i][k]+taus[j+1][i][k]+taus[j+1][i][k+1]+taus[j][i][k+1]);
tausipkp[j][i][k]=0.25*(taus[j][i][k]+taus[j][i+1][k]+taus[j][i+1][k+1]+taus[j][i][k+1]);
mod_av->tausipjp[j][i][k]=0.25*(mod->taus[j][i][k]+mod->taus[j][i+1][k]+mod->taus[j+1][i+1][k]+mod->taus[j+1][i][k]);
mod_av->tausjpkp[j][i][k]=0.25*(mod->taus[j][i][k]+mod->taus[j+1][i][k]+mod->taus[j+1][i][k+1]+mod->taus[j][i][k+1]);
mod_av->tausipkp[j][i][k]=0.25*(mod->taus[j][i][k]+mod->taus[j][i+1][k]+mod->taus[j][i+1][k+1]+mod->taus[j][i][k+1]);
}
rjp[j][i][k]=0.5*(rho[j][i][k]+rho[j+1][i][k]);
rkp[j][i][k]=0.5*(rho[j][i][k]+rho[j][i][k+1]);
rip[j][i][k]=0.5*(rho[j][i][k]+rho[j][i+1][k]);
mod_av->rjp[j][i][k]=0.5*(mod->rho[j][i][k]+mod->rho[j+1][i][k]);
mod_av->rkp[j][i][k]=0.5*(mod->rho[j][i][k]+mod->rho[j][i][k+1]);
mod_av->rip[j][i][k]=0.5*(mod->rho[j][i][k]+mod->rho[j][i+1][k]);
}
......
......@@ -37,8 +37,7 @@ void err2(char errformat[],char errfilename[]){
err(outtxt);
}
void checkfd(FILE *fp, float *** prho, float *** ppi, float *** pu,
float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **recpos, int ntr){
void checkfd(FILE *fp, st_model *mod, float **srcpos, int nsrc, int **recpos, int ntr){
extern float DX, DY, DZ, DT, TS, TIME, TSNAP2;
extern float XREC1, XREC2, YREC1, YREC2, ZREC1, ZREC2;
......@@ -80,12 +79,12 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
for (k=1+nfw;k<=(nz-nfw);k++){
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else{
if(L){
c=sqrt(pu[j][i][k]*(1.0+L*ptaus[j][i][k])/prho[j][i][k]);}
else c=sqrt(pu[j][i][k]/prho[j][i][k]);
c=sqrt(mod->u[j][i][k]*(1.0+L*mod->taus[j][i][k])/mod->rho[j][i][k]);}
else c=sqrt(mod->u[j][i][k]/mod->rho[j][i][k]);
}
if (cmax_s<c) cmax_s=c;
......@@ -93,13 +92,13 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
frequency of the source */
sum=0.0;
for (l=1;l<=L;l++){
ts=DT/peta[l];
sum=sum+((w*w*ptaus[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
ts=DT/mod->eta[l];
sum=sum+((w*w*mod->taus[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
}
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else
c=sqrt(pu[j][i][k]*(1.0+sum)/prho[j][i][k]);
c=sqrt(mod->u[j][i][k]*(1.0+sum)/mod->rho[j][i][k]);
if ((c>cwater)&&(cmin_s>c)) cmin_s=c;
}
}
......@@ -112,24 +111,24 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
for (k=1+nfw;k<=(nz-nfw);k++){
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else{
if(L) c=sqrt(ppi[j][i][k]*(1.0+L*ptaup[j][i][k])/prho[j][i][k]);
else c=sqrt(ppi[j][i][k]/prho[j][i][k]);
if(L) c=sqrt(mod->pi[j][i][k]*(1.0+L*mod->taup[j][i][k])/mod->rho[j][i][k]);
else c=sqrt(mod->pi[j][i][k]/mod->rho[j][i][k]);
}
if (cmax_p<c) cmax_p=c;
/* find minimum model phase velocity of P-waves at center
frequency of the source */
sum=0.0;
for (l=1;l<=L;l++){
ts=DT/peta[l];
sum=sum+((w*w*ptaup[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
ts=DT/mod->eta[l];
sum=sum+((w*w*mod->taup[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
}
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else
c=sqrt(ppi[j][i][k]*(1.0+sum)/prho[j][i][k]);
c=sqrt(mod->pi[j][i][k]*(1.0+sum)/mod->rho[j][i][k]);
if (cmin_p>c) cmin_p=c;
}
}
......
......@@ -26,10 +26,7 @@
#include "fd.h"
void comm_ini(float *** bufferlef_to_rig, float *** bufferrig_to_lef,
float *** buffertop_to_bot, float *** bufferbot_to_top,
float *** bufferfro_to_bac, float *** bufferbac_to_fro,
MPI_Request *req_send, MPI_Request *req_rec){
void comm_ini(st_buffer *buffer, MPI_Request *req_send, MPI_Request *req_rec){
extern int NX, NY, NZ, INDEX[7], FDORDER;
......@@ -50,22 +47,24 @@ MPI_Request *req_send, MPI_Request *req_rec){
nf2=nf1-1;
MPI_Bsend_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&bufferfro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&bufferbac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffertop_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&bufferbot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
MPI_Bsend_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
/* initialising of receive of buffer arrays. Same arrays for send and receive
are used. Thus, before starting receive, it is necessary to check if
Bsend has copied data into local buffers, i.e. has completed.
*/
MPI_Recv_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&bufferfro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&bufferbac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffertop_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&bufferbot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
MPI_Recv_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
}
......@@ -7,10 +7,7 @@
#include "fd.h"
void comm_ini_s(float *** bufferlef_to_rig, float *** bufferrig_to_lef,
float *** buffertop_to_bot, float *** bufferbot_to_top,
float *** bufferfro_to_bac, float *** bufferbac_to_fro,
MPI_Request *req_send, MPI_Request *req_rec){
void comm_ini_s(st_buffer *buffer, MPI_Request *req_send, MPI_Request *req_rec){
extern int NX, NY, NZ, INDEX[7], FDORDER;
......@@ -31,22 +28,22 @@ MPI_Request *req_send, MPI_Request *req_rec){
nf2=nf1-1;
MPI_Bsend_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&bufferfro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&bufferbac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffertop_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&bufferbot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
MPI_Bsend_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
/* initialising of receive of buffer arrays. Same arrays for send and receive
are used. Thus, before starting receive, it is necessary to check if
Bsend has copied data into local buffers, i.e. has completed.
*/
MPI_Recv_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&bufferfro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&bufferbac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffertop_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&bufferbot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
MPI_Recv_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
}
......@@ -24,7 +24,7 @@
#include "fd.h"
void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***grad3, float ***gradprior1, float ***gradprior2, float ***gradprior3, float ***gradprior4, float ***gradprior5, float ***gradprior6,float *beta, int iteration,int cdf){
void conjugate(int nx,int ny,int nz, st_gradient *grad, st_gradient *grad_prior1,st_gradient *grad_prior2, float *beta, int iteration,int cdf){
extern FILE *FP;
......@@ -43,12 +43,12 @@ void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***gra
for (j=1;j<=ny;j++){
for (i=1;i<=nx;i++){
for (k=1;k<=nz;k++){
gradprior1[j][i][k]=grad1[j][i][k];
gradprior2[j][i][k]=grad2[j][i][k];
gradprior3[j][i][k]=grad3[j][i][k];
gradprior4[j][i][k]=grad1[j][i][k];
gradprior5[j][i][k]=grad2[j][i][k];
gradprior6[j][i][k]=grad3[j][i][k];
grad_prior1->vp[j][i][k]=grad->vp[j][i][k];
grad_prior1->vs[j][i][k]=grad->vs[j][i][k];
grad_prior1->rho[j][i][k]=grad->rho[j][i][k];
grad_prior2->vp[j][i][k]=grad->vp[j][i][k];
grad_prior2->vs[j][i][k]=grad->vs[j][i][k];
grad_prior2->rho[j][i][k]=grad->rho[j][i][k];
}
}
}
......@@ -59,17 +59,17 @@ void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***gra
for (j=1;j<=ny;j++){
for (i=1;i<=nx;i++){
for (k=1;k<=nz;k++){
betadum[0]+=grad1[j][i][k]*(grad1[j][i][k]-gradprior1[j][i][k]);
betadum[1]+=gradprior1[j][i][k]*gradprior1[j][i][k];
betadum[2]+=grad2[j][i][k]*(grad2[j][i][k]-gradprior2[j][i][k]);
betadum[3]+=gradprior2[j][i][k]*gradprior2[j][i][k];
betadum[4]+=grad3[j][i][k]*(grad3[j][i][k]-gradprior3[j][i][k]);
betadum[5]+=gradprior3[j][i][k]*gradprior3[j][i][k];
betadum[0]+=grad->vp[j][i][k]*(grad->vp[j][i][k]-grad_prior1->vp[j][i][k]);
betadum[1]+=grad_prior1->vp[j][i][k]*grad_prior1->vp[j][i][k];
betadum[2]+=grad->vs[j][i][k]*(grad->vs[j][i][k]-grad_prior1->vs[j][i][k]);
betadum[3]+=grad_prior1->vs[j][i][k]*grad_prior1->vs[j][i][k];
betadum[4]+=grad->rho[j][i][k]*(grad->rho[j][i][k]-grad_prior1->rho[j][i][k]);
betadum[5]+=grad_prior1->rho[j][i][k]*grad_prior1->rho[j][i][k];
/*save preconditioned gradient for next iteration*/
gradprior1[j][i][k]=grad1[j][i][k];
gradprior2[j][i][k]=grad2[j][i][k];
gradprior3[j][i][k]=grad3[j][i][k];
grad_prior1->vp[j][i][k]=grad->vp[j][i][k];
grad_prior1->vs[j][i][k]=grad->vs[j][i][k];
grad_prior1->rho[j][i][k]=grad->rho[j][i][k];
}
}
}
......@@ -94,14 +94,14 @@ void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***gra
for (j=1;j<=ny;j++){
for (i=1;i<=nx;i++){
for (k=1;k<=nz;k++){
grad1[j][i][k]=gradprior4[j][i][k]*beta[0]+grad1[j][i][k];
grad2[j][i][k]=gradprior5[j][i][k]*beta[1]+grad2[j][i][k];
grad3[j][i][k]=gradprior6[j][i][k]*beta[2]+grad3[j][i][k];
grad->vp[j][i][k]=grad_prior2->vp[j][i][k]*beta[0]+grad->vp[j][i][k];
grad->vs[j][i][k]=grad_prior2->vs[j][i][k]*beta[1]+grad->vs[j][i][k];
grad->rho[j][i][k]=grad_prior2->rho[j][i][k]*beta[2]+grad->rho[j][i][k];
/*save conjugate gradient for next iteration*/
gradprior4[j][i][k]=grad1[j][i][k];
gradprior5[j][i][k]=grad2[j][i][k];
gradprior6[j][i][k]=grad3[j][i][k];
grad_prior2->vp[j][i][k]=grad->vp[j][i][k];
grad_prior2->vs[j][i][k]=grad->vs[j][i][k];
grad_prior2->rho[j][i][k]=grad->rho[j][i][k];
}
}
}
......
This diff is collapsed.
......@@ -2,46 +2,60 @@
* Copyright (C) 2015 For the list of authors, see file AUTHORS.
*
* This file is part of IFOS3D.
*
*
* IFOS3D is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
*
* IFOS3D is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
*
* You should have received a copy of the GNU General Public License
* along with IFOS3D. See file COPYING and/or
* along with IFOS3D. See file COPYING and/or
* <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
Copying of model parmeters into a testmodel
Copying of model parmeters into a testmodel
S. Butzer 2013
* ----------------------------------------------------------------------*/
#include "fd.h"
void cpmodel(int nx, int ny, int nz, float ***rho, float ***pi, float ***u,float *** testrho, float *** testpi, float *** testu){
int j,i,k;
for (j=1;j<=ny;j++){
for (i=1;i<=nx;i++){
for (k=1;k<=nz;k++){
testrho[j][i][k]=0.0;
testrho[j][i][k]=rho[j][i][k];
testu[j][i][k]=0.0;
testu[j][i][k]=u[j][i][k];
testpi[j][i][k]=0.0;
testpi[j][i][k]=pi[j][i][k];
void cpmodel(int nx, int ny, int nz, st_model *mod,st_model *testmod) {
extern int L;
int j,i,k,l;
for (j=1; j<=ny; j++) {
for (i=1; i<=nx; i++) {
for (k=1; k<=nz; k++) {
testmod->rho[j][i][k]=0.0;
testmod->rho[j][i][k]=mod->rho[j][i][k];
testmod->u[j][i][k]=0.0;
testmod->u[j][i][k]=mod->u[j][i][k];
testmod->pi[j][i][k]=0.0;
testmod->pi[j][i][k]=mod->pi[j][i][k];
if (L) {
testmod->taup[j][i][k]=0.0;
testmod->taup[j][i][k]=mod->taup[j][i][k];
testmod->taus[j][i][k]=0.0;
testmod->taus[j][i][k]=mod->taus[j][i][k];
}
}
}
}
\ No newline at end of file
}
if (L) {
for (l=1; l<=L; l++) {
testmod->eta[l]=0.0;
testmod->eta[l]=mod->eta[l];
}
}
}
......@@ -25,7 +25,7 @@
#include "fd.h"
void discfourier(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, int nt, float *** vx, float *** vy, float *** vz, float **** Fvx, float **** Fvy,float **** Fvz, float **** Fvxi, float **** Fvyi,float **** Fvzi,float *finv, int nf, int ntast, int pshot1, int back){
void discfourier(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, int nt, st_velocity *vel, st_freq_velocity *fourier_vel,float *finv, int nf, int ntast, int pshot1, int back){
extern float DT,TIME;
......@@ -47,13 +47,13 @@ void discfourier(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, int nt, f
for (i=nx1;i<=nx2;i++){
for (k=nz1;k<=nz2;k++){
Fvx[m][j][i][k]+=vx[j][i][k]*trig1;
Fvy[m][j][i][k]+=vy[j][i][k]*trig1;
Fvz[m][j][i][k]+=vz[j][i][k]*trig1;
fourier_vel->Fvx_re[m][j][i][k]+=vel->vx[j][i][k]*trig1;
fourier_vel->Fvy_re[m][j][i][k]+=vel->vy[j][i][k]*trig1;
fourier_vel->Fvz_re[m][j][i][k]+=vel->vz[j][i][k]*trig1;
Fvxi[m][j][i][k]+=vx[j][i][k]*trig2;
Fvyi[m][j][i][k]+=vy[j][i][k]*trig2;
Fvzi[m][j][i][k]+=vz[j][i][k]*trig2;
fourier_vel->Fvx_im[m][j][i][k]+=vel->vx[j][i][k]*trig2;
fourier_vel->Fvy_im[m][j][i][k]+=vel->vy[j][i][k]*trig2;
fourier_vel->Fvz_im[m][j][i][k]+=vel->vz[j][i][k]*trig2;
}
}
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -25,7 +25,7 @@
#include "fd.h"
void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***grad1, float ***grad2,float ***grad3, float ***hess1, float ***hess2, float ***hess3, float finv, int iteration){
void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, st_gradient *grad, st_hessian *hessian, float finv, int iteration){
extern float WATER_HESS[3], VP0, VS0, RHO0;
extern int NXG, NYG, NZG, FW, NX, NY, NZ;
......@@ -44,12 +44,12 @@ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***g
for (j=1;j<=ny2;j++){
for (i=1;i<=nx2;i++){
for (k=1;k<=nz2;k++){
grad1[j][i][k]=grad1[j][i][k]*VP0;
grad2[j][i][k]=grad2[j][i][k]*VS0;
grad3[j][i][k]=grad3[j][i][k]*RHO0;
/*hess1[j][i][k]=hess1[j][i][k]*pow(vp0,2.0);
hess2[j][i][k]=hess2[j][i][k]*pow(vs0,2.0);
hess3[j][i][k]=hess3[j][i][k]*pow(rho0,2.0);*/
grad->vp[j][i][k]=grad->vp[j][i][k]*VP0;
grad->vs[j][i][k]=grad->vs[j][i][k]*VS0;
grad->rho[j][i][k]=grad->rho[j][i][k]*RHO0;
/*hessian->vp[j][i][k]=hessian->vp[j][i][k]*pow(vp0,2.0);
hessian->vs[j][i][k]=hessian->vs[j][i][k]*pow(vs0,2.0);
hessian->rho[j][i][k]=hessian->rho[j][i][k]*pow(rho0,2.0);*/
}
}
}
......@@ -59,9 +59,9 @@ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***g
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
for (k=nz1;k<=nz2;k++){
wl[0]+=log10(hess1[j][i][k]);
wl[1]+=log10(hess2[j][i][k]);
wl[2]+=log10(hess3[j][i][k]);
wl[0]+=log10(hessian->vp[j][i][k]);
wl[1]+=log10(hessian->vs[j][i][k]);
wl[2]+=log10(hessian->rho[j][i][k]);
}
}
}
......@@ -79,9 +79,9 @@ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***g
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
grad1[j][i][k]=grad1[j][i][k]/(hess1[j][i][k]+WATER_HESS[0]);
grad2[j][i][k]=grad2[j][i][k]/(hess2[j][i][k]+WATER_HESS[1]);
grad3[j][i][k]=grad3[j][i][k]/(hess3[j][i][k]+WATER_HESS[2]);
grad->vp[j][i][k]=grad->vp[j][i][k]/(hessian->vp[j][i][k]+WATER_HESS[0]);
grad->vs[j][i][k]=grad->vs[j][i][k]/(hessian->vs[j][i][k]+WATER_HESS[1]);
grad->rho[j][i][k]=grad->rho[j][i][k]/(hessian->rho[j][i][k]+WATER_HESS[2]);
}
}
}
......
This diff is collapsed.
......@@ -23,7 +23,7 @@
--------------------------------------------------------------------------*/
#include "fd.h"
void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, float **bfgsmod, float **bfgsgrad, int iteration){
void lbfgs(st_gradient *grad, float *bfgsscale, float **bfgsmod, float **bfgsgrad, int iteration){
int m=0,v=0,w=0;
int i,j,k,l;
......@@ -53,8 +53,8 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
l++;
bfgsgrad[w][l]+=-grad1[j][i][k];
q[l]=grad1[j][i][k];
bfgsgrad[w][l]+=-grad->vp[j][i][k];
q[l]=grad->vp[j][i][k];
}
}
}
......@@ -64,8 +64,8 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
l++;
bfgsgrad[w][l]+=-grad2[j][i][k];
q[l]=grad2[j][i][k];
bfgsgrad[w][l]+=-grad->vs[j][i][k];
q[l]=grad->vs[j][i][k];
}
}
}
......@@ -159,8 +159,8 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
l++; /*w=(j-1)*Nx*Nz+(i-1)*Nz+k*/
bfgsgrad[w][l]=grad1[j][i][k];
grad1[j][i][k]=q[l];
bfgsgrad[w][l]=grad->vp[j][i][k];
grad->vp[j][i][k]=q[l];
}
}
}
......@@ -170,8 +170,8 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
l++; /*w=(j-1)*Nx*Nz+(i-1)*Nz+k*/
bfgsgrad[w][l]=grad2[j][i][k];
grad2[j][i][k]=q[l];
bfgsgrad[w][l]=grad->vs[j][i][k];
grad->vs[j][i][k]=q[l];
}
}
}
......
......@@ -22,7 +22,7 @@
--------------------------------------------------------------------------------------*/
#include "fd.h"
void lbfgs_savegrad(float ***grad1, float ***grad2, float ***grad3,float **bfgsgrad1){
void lbfgs_savegrad(st_gradient *grad,float **bfgsgrad1){
int d,i,j,k;
extern int NX,NY, NZ, NUMPAR;
......@@ -33,8 +33,8 @@ void lbfgs_savegrad(float ***grad1, float ***grad2, float ***grad3,float **bfgsg
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
d++;
bfgsgrad1[1][d]=+grad1[j][i][k];
/*bfgsgrad3[1][d]=+grad3[j][i][k];*/
bfgsgrad1[1][d]=+grad->vp[j][i][k];
/*bfgsgrad3[1][d]=+grad->rho[j][i][k];*/
}
}
......@@ -45,7 +45,7 @@ void lbfgs_savegrad(float ***grad1, float ***grad2, float ***grad3,float **bfgsg
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
d++;
bfgsgrad1[1][d]=+grad2[j][i][k];
bfgsgrad1[1][d]=+grad->vs[j][i][k];
}
}
}
......
......@@ -28,8 +28,7 @@
#include "fd.h"
void matcopy(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup){
void matcopy(st_model *mod){
extern int MYID, NX, NY, NZ, INDEX[7],L,VERBOSE;
......@@ -67,13 +66,13 @@ float *** taus, float *** taup){
for (k=0;k<=NZ+1;k++){
/* storage of top of local volume into buffer */
buffertop_to_bot[i][k][1] = rho[1][i][k];
buffertop_to_bot[i][k][2] = pi[1][i][k];
buffertop_to_bot[i][k][3] = u[1][i][k];
buffertop_to_bot[i][k][1] = mod->rho[1][i][k];
buffertop_to_bot[i][k][2] = mod->pi[1][i][k];
buffertop_to_bot[i][k][3] = mod->u[1][i][k];
if(L){
buffertop_to_bot[i][k][4] = taus[1][i][k];
buffertop_to_bot[i][k][5] = taup[1][i][k];
buffertop_to_bot[i][k][4] = mod->taus[1][i][k];
buffertop_to_bot[i][k][5] = mod->taup[1][i][k];
}
}
}
......@@ -84,13 +83,13 @@ float *** taus, float *** taup){
/* storage of bottom of local volume into buffer */
bufferbot_to_top[i][k][1] = rho[NY][i][k];
bufferbot_to_top[i][k][2] = pi[NY][i][k];
bufferbot_to_top[i][k][3] = u[NY][i][k];
bufferbot_to_top[i][k][1] = mod->rho[NY][i][k];
bufferbot_to_top[i][k][2] = mod->pi[NY][i][k];
bufferbot_to_top[i][k][3] = mod->u[NY][i][k];
if(L){
bufferbot_to_top[i][k][4] = taus[NY][i][k];
bufferbot_to_top[i][k][5] = taup[NY][i][k];
bufferbot_to_top[i][k][4] = mod->taus[NY][i][k];
bufferbot_to_top[i][k][5] = mod->taup[NY][i][k];
}
}
......@@ -111,12 +110,12 @@ float *** taus, float *** taup){
for (i=0;i<=NX+1;i++){
for (k=0;k<=NZ+1;k++){
rho[NY+1][i][k] = buffertop_to_bot[i][k][1];
pi[NY+1][i][k] = buffertop_to_bot[i][k][2];
u[NY+1][i][k] = buffertop_to_bot[i][k][3];
mod->rho[NY+1][i][k] = buffertop_to_bot[i][k][1];
mod->pi[NY+1][i][k] = buffertop_to_bot[i][k][2];
mod->u[NY+1][i][k] = buffertop_to_bot[i][k][3];
if(L){
taus[NY+1][i][k] = buffertop_to_bot[i][k][4];
taup[NY+1][i][k] = buffertop_to_bot[i][k][5];
mod->taus[NY+1][i][k] = buffertop_to_bot[i][k][4];
mod->taup[NY+1][i][k] = buffertop_to_bot[i][k][5];
}
......@@ -126,12 +125,12 @@ float *** taus, float *** taup){
for (i=0;i<=NX+1;i++){
for (k=0;k<=NZ+1;k++){
rho[0][i][