Commit 3363b15f authored by Florian Wittkamp's avatar Florian Wittkamp

Small clean up

-Remove src_relicts
-Clean up of the Makefile in src/
parent 52828193
172.22.135.9 slots=8
172.22.135.12 slots=8
172.22.135.15 slots=8
......@@ -12,31 +12,26 @@ MODEL_AC = ../genmod/1D_linear_gradient_ac.c
MODEL_VAC = ../genmod/1D_linear_gradient_viscac.c
EXEC= ../bin
# Compiler (LAM: CC=hcc, CRAY T3E: CC=cc)
# ON Linux cluster running LAM
#CC=hcc
#LFLAGS=-lm -lmpi
#CFLAGS=-Wall -O4
# On CRAY T3E
# CC=cc
# On CHIC
CC=mpicc
# Description:
# CC = Compiler
# LFLAGS = Linker flag
# CFLAGS = Compiler flag
# LINUX with OpenMPI / IntelMPI and INTEL Compiler
# Use icc whenever possible, this will be much faster than gcc
CC=mpiicc
LFLAGS=-lm -lcseife -lstfinv -laff -lfourierxx -lfftw3 -lstdc++
CFLAGS=-O3
SFLAGS=-L./../contrib/libcseife -L./../contrib/bin
IFLAGS=-I./../contrib/libcseife -I./../contrib/header -I.
# On HLRN system
#CC=mpcc
#LFLAGS=-lm
# LINUX with OpenMPI / IntelMPI and GCC Compiler
#CC=mpicc
#LFLAGS=-lm -lcseife -lstfinv -laff -lfourierxx -lfftw3 -lstdc++
#CFLAGS=-O3
#SFLAGS=-L./../contrib/libcseife -L./../contrib/bin
#IFLAGS=-I./../contrib/libcseife -I./../contrib/header -I.
# ALTIX
#CC=icc
#CFLAGS=-mp -O3 -ipo
#LFLAGS=-lmpi -lm -i-static
# after this line, no further editing should be necessary
# --------------------------------------------------------
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* solve linear matrix equation using LU decomposition
*
*
* Daniel Koehn
* last update 9.11.2007
*
* ---------------------------------------------------------------------*/
#include "fd.h"
double LU_decomp(double **A, double *x, double *b,int n){
int i,j,k,s;
double suma, *y;
y = dvector(1,n);
for (j=1;j<=n;j++){
for (i=1;i<=j;i++){
if(i==1){
A[i][j] = A[i][j];
}else{
suma = 0.0;
for (k=1;k<=i-1;k++){
suma = suma + A[i][k]*A[k][j];
}
A[i][j] = A[i][j] - suma;
}
}
if (j<n){
for (s=1;s<=(n-j);s++){
i = j + s;
if (j==1){
A[i][j] = A[i][j]/A[j][j];
}else{
suma = 0.0;
for (k=1;k<=(j-1);k++){
suma = suma + A[i][k]*A[k][j];
}
A[i][j] = (A[i][j] - suma)/A[j][j];
}
}
}
}
y[1] = b[1];
for (i=2;i<=n;i++){
suma = 0.0;
for (j=1;j<=(i-1);j++){
suma = suma + A[i][j]*y[j];
}
y[i] = b[i] - suma;
}
i = n;
j = n;
x[i] = y[i]/A[i][j];
for (s=1;s<=(n-1);s++){
i = n - s;
suma = 0.0;
for (j=i+1;j<=n;j++){
suma = suma + A[i][j]*x[j];
}
x[i] = (y[i] - suma)/A[i][i];
}
free_dvector(y,1,n);
return 0;
}
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*
* Create dissipative boundarie around the model grid
* The dissipative coefficients are stored in the 2-D array
* absorb_coeff. The interior of the model is weighted by the
* coefficient 1. In the absorbing frame the coefficients
* are less than one. Coefficients are computed using
* exponential damping (see Cerjan et al., 1985,
* Geophysics, 50, 705-708)
*
* last update K. Witte 22.12.01 :-)
*/
#include "fd.h"
void absorb(float ** absorb_coeff)
{
/* extern variables */
extern float DH, FW, DAMPING;
extern int FREE_SURF, NX, NY, BOUNDARY;
extern int NPROCX, NPROCY, MYID, POS[3];
extern FILE *FP;
/* local variables */
int i, j, ifw, ii, jj, xb, yb, xe, ye;
float amp, a, *coeff;
/*char modfile[STRING_SIZE];*/
if (MYID==0)
{
fprintf(FP,"\n **Message from absorb (printed by PE %d):\n",MYID);
fprintf(FP," Coefficcients for absorbing frame are now calculated.\n");
fprintf(FP," Width of dissipative frame (meter)= %f\n",FW);
fprintf(FP," Percentage of exponential damping = %5.2f\n",DAMPING);
}
amp=1.0-DAMPING/100.0; /* amplitude at the edge of the numerical grid */
ifw=iround(FW/DH); /* frame width in gridpoints */
coeff=vector(1,ifw);
a=sqrt(-log(amp)/((ifw-1)*(ifw-1)));
for (i=1;i<=ifw;i++)
coeff[i]=exp(-(a*a*(ifw-i)*(ifw-i)));
if (MYID==0)
{
fprintf(FP," Table of coefficients \n # \t coeff \n");
/*printf(" ifw=%d \t a=%f amp=%f \n", ifw,a,amp); */
for (i=1;i<=ifw;i++)
fprintf(FP," %d \t %5.3f \n", i, coeff[i]);
}
/* initialize array of coefficients with one */
for (j=1;j<=NY;j++)
for (i=1;i<=NX;i++) absorb_coeff[j][i]=1.0;
/* compute coefficients for left and right grid boundaries (x-direction) */
if ((!BOUNDARY) && (POS[1]==0))
{
yb=1; ye=NY;
for (i=1;i<=ifw;i++)
{
if ((POS[2]==0) && (!(FREE_SURF))) yb=i;
if (POS[2]==NPROCY-1) ye=NY-i+1;
for (j=yb;j<=ye;j++)
absorb_coeff[j][i]=coeff[i];
}
}
if ((!BOUNDARY) && (POS[1]==NPROCX-1))
{
yb=1; ye=NY;
for (i=1;i<=ifw;i++){
ii=NX-i+1;
if ((POS[2]==0) && (!(FREE_SURF))) yb=i;
if (POS[2]==NPROCY-1) ye=NY-i+1;
for (j=yb;j<=ye;j++)
absorb_coeff[j][ii]=coeff[i];
}
}
/* compute coefficients for top and bottom grid boundaries (y-direction) */
if ((POS[2]==0) && (!(FREE_SURF)))
{
xb=1; xe=NX;
for (j=1;j<=ifw;j++)
{
if ((!BOUNDARY) && (POS[1]==0)) xb=j;
if ((!BOUNDARY) && (POS[1]==NPROCX-1)) xe=NX-j+1;
for (i=xb;i<=xe;i++)
absorb_coeff[j][i]=coeff[j];
}
}
if (POS[2]==NPROCY-1)
{
xb=1; xe=NX;
for (j=1;j<=ifw;j++)
{
jj=NY-j+1;
if ((!BOUNDARY) && (POS[1]==0)) xb=j;
if ((!BOUNDARY) && (POS[1]==NPROCX-1)) xe=NX-j+1;
for (i=xb;i<=xe;i++)
absorb_coeff[jj][i]=coeff[j];
}
}
/* sprintf(modfile,"absorb_coeff.bin");
writemod(modfile,absorb_coeff,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
*/
free_vector(coeff,1,ifw);
}
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* calculate step length for material parameter update
*
* Daniel Koehn
* last update 9.11.2007
*
* ---------------------------------------------------------------------*/
#include "fd.h"
void calc_mat_change(float ** waveconv, float ** waveconv_rho, float ** waveconv_u, float ** rho, float ** rhonp1, float ** pi, float ** pinp1, float ** u, float ** unp1, int iter,
int epstest, int INVMAT, float eps_scale_vp, float eps_scale_vs){
/*--------------------------------------------------------------------------*/
FILE *FP1;
/* extern variables */
extern float DH, DT;
extern float EPSILON, MUN, EPSILON_u;
extern int NX, NY, NXG, NYG, POS[3], MYID;
/* local variables */
float Rho, Vp, Vs, Vsnp1, Vpnp1, x, y, undf, r, pi0, K, mu;
float dpi, pimax, rhomax, umax,gradmax, gradmax_rho, gradmax_u, epsilon1, pimaxr, gradmaxr, umaxr, gradmaxr_u;
float VPUPPERLIM, VSLOWERLIM;
int i, j, ii, jj;
char modfile[STRING_SIZE];
VPUPPERLIM = 2200.0;
VSLOWERLIM = 800.0;
/* find maximum of pi and gradient waveconv */
if((INVMAT==1)||(INVMAT==0)){
pimax = 0.0;
gradmax = 0.0;
for (i=1;i<=NX;i++){
for (j=1;j<=NY;j++){
Vp = sqrt((pi[j][i] + 2.0 * u[j][i])/rho[j][i]);
if(Vp>pimax){pimax=Vp;}
if((i*j == 1) || (gradmax == 0.0)) {
gradmax = fabs(waveconv[j][i]);
} else {
if(fabs(waveconv[j][i]) > gradmax){
gradmax = fabs(waveconv[j][i]);
}
}
}}}
/* find maximum of u and gradient waveconv_u */
if((INVMAT==3)||(INVMAT==0)){
umax = 0.0;
gradmax_u = 0.0;
for (i=1;i<=NX;i++){
for (j=1;j<=NY;j++){
Vs = sqrt(u[j][i]/rho[j][i]);
if(Vs>umax){umax=Vs;}
if((i*j == 1) || (gradmax_u == 0.0)) {
gradmax_u = fabs(waveconv_u[j][i]);
} else {
if(fabs(waveconv_u[j][i]) > gradmax_u){
gradmax_u = fabs(waveconv_u[j][i]);
}
}
}}}
/* calculate test step length epsilon for pi */
if(INVMAT==0){
/*eps_scale_vp = 150.0;*/
pimaxr=0.0;
gradmaxr=0.0;
MPI_Allreduce(&pimax, &pimaxr, 1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&gradmax,&gradmaxr,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
EPSILON = (pimaxr/(eps_scale_vp*(gradmaxr)));
epsilon1 = EPSILON;
MPI_Allreduce(&EPSILON,&epsilon1,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
if (MYID==0) EPSILON=epsilon1;
exchange_par();
}
if(INVMAT==0){
umaxr=0.0;
gradmaxr_u=0.0;
MPI_Allreduce(&umax,&umaxr, 1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&gradmax_u,&gradmaxr_u,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
EPSILON_u = (umaxr/(eps_scale_vp*gradmaxr_u));
epsilon1 = EPSILON_u;
MPI_Allreduce(&EPSILON_u,&epsilon1,1,MPI_FLOAT,MPI_MIN,MPI_COMM_WORLD);
if (MYID==0) EPSILON_u=epsilon1;
exchange_par();
}
if(epstest==1){
/* loop over local grid */
for (i=1;i<=NX;i++){
for (j=1;j<=NY;j++){
/* update material parameters */
if(INVMAT==0){
Vs = sqrt(u[j][i]/rho[j][i]);
Vp = sqrt((pi[j][i] + 2.0 * u[j][i])/rho[j][i]);
Vp = Vp + EPSILON*waveconv[j][i];
Vs = Vs + EPSILON_u*waveconv_u[j][i];
/*Vs = Vp/sqrt(3.0);*/
if (Vp>VPUPPERLIM) Vp = VPUPPERLIM;
if (Vs<VSLOWERLIM) Vs = VSLOWERLIM;
/*rho[j][i] = 1000.0*0.31*pow(Vp,(1.0/4.0));*/
pi[j][i] = Vp*Vp*rho[j][i] - (2.0 * Vs*Vs*rho[j][i]);
u[j][i] = Vs*Vs*rho[j][i];
}
/*if((INVMAT==2)){
rho[j][i] = rho[j][i] - EPSILON_rho*waveconv_rho[j][i];}*/
/*if(INVMAT==3){
mu = u[j][i];
mu = mu + EPSILON_u*waveconv_u[j][i];
u[j][i] = mu;}*/
/* output of Vp and Vs */
unp1[j][i] = sqrt(u[j][i]/rho[j][i]);
pinp1[j][i] = sqrt((pi[j][i]+2.0*u[j][i])/rho[j][i]);
}
}
sprintf(modfile,"model/waveform_test_model_vp_it_%d.bin",iter);
writemod(modfile,pinp1,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
sprintf(modfile,"model/waveform_test_model_vs_it_%d.bin",iter);
writemod(modfile,unp1,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
sprintf(modfile,"model/waveform_test_model_rho_it_%d.bin",iter);
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
}
}
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* calculate step length for material parameter update
*
* waveconv = conjugated gradient direction
* gradp = preconditionated gradient
*
* Daniel Koehn
* last update 9.11.2007
*
* ---------------------------------------------------------------------*/
#include "fd.h"
float calc_opt_step_test(float * L2t, float ** waveconv, float ** gradg, float * epst, int sws, float C_vp){
extern int NX, NY, IDX, IDY, MYID;
extern float EPSILON, EPSILON_u, EPSILON_rho;
int i, j, n;
float opteps, H1, H1sum, H2, H2sum, etas, H1max, H1maxsum;
/* calculate optimal step size after Tarantola (1986)*/
/*H1max=0.0;
H1=0.0;
H2=0.0;
for (i=1;i<=NX;i=i+IDX){
for (j=1;j<=NY;j=j+IDY){
H1 += 1e+39 * waveconv[j][i]*gradg[j][i];
H2 += waveconv[j][i]*(1.0/C_vp)*waveconv[j][i];
if(fabs(H1)>H1max){H1max=H1;}
}
}
H1sum=0.0;
MPI_Allreduce(&H1,&H1sum,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
H2sum=0.0;
MPI_Allreduce(&H2,&H2sum,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
H1maxsum=0.0;
MPI_Allreduce(&H1max,&H1maxsum,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);*/
etas = epst[sws];
/* calculate optimal step length for Vp update */
opteps = etas * L2t[sws];
if((opteps/etas) > 2.0){opteps = 2.0 * etas;}
printf("MYID = %d \t L2t[sws] = %e \t epst[sws] = %e \t sws = %d \t opteps = %e \t C_vp = %e \n ",MYID,L2t[sws],epst[sws],sws,opteps,C_vp);
return opteps;
}
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Write seismograms to disk
* last update 19/01/02, T. Bohlen
* ----------------------------------------------------------------------*/
#include "fd.h"
#include "segy.h"
void conv_FD(float * temp_TS, float * temp_TS1, float * temp_conv, int ns){
/* declaration of extern variables */
extern int NDT, MYID, TAPERLENGTH;
extern char DATA_DIR[STRING_SIZE];
extern float TIME, DH, DT, REFREC[4];
const float xshift=800.0, yshift=800.0;
/* declaration of local variables */
int i,j, h, nfreq, fa, fb, npad, itr;
int tracl1 ;
float xr, yr, x, y, dump, a, times;
float damping, amp, maximum;
float *window, *amp1;
float *trace_real, *trace_complex, *trace_real1, *trace_complex1;
float *trace_real_conv, *trace_complex_conv;
double npadd;
/* declaration of local variables */
int if1, if2, if3, if4;
npad = (int)(pow(2.0, ceil(log((double)(ns))/log(2.0))+2.0) ); /* ns -> npad for usage in FFT*/
npadd = (double)npad;
trace_real = vector(1,npad);
trace_complex = vector(1,npad);
trace_real1 = vector(1,npad);
trace_complex1 = vector(1,npad);
trace_real_conv = vector(1,npad);
trace_complex_conv = vector(1,npad);
window = vector(1,(npad/2));
/* calculate TIMEs after application of zero pad*/
times=npad*(TIME/ns);
for(j=1;j<=ns;j++){
trace_real[j] = temp_TS[j];
trace_complex[j] = 0.0;
trace_real1[j] = temp_TS1[j];
trace_complex1[j] = 0.0;
}
/* apply fft */
FFT(1,(unsigned long)(log(npadd)/log(2.0)),trace_real,trace_complex);
FFT(1,(unsigned long)(log(npadd)/log(2.0)),trace_real1,trace_complex1);
/* multiplication of complex vectors */
for(i=1;i<=npad;i++){
trace_real_conv[i] = (trace_real[i] * trace_real1[i]) - (trace_complex[i] * trace_complex1[i]);
trace_complex_conv[i] = (trace_real[i] * trace_complex1[i]) + (trace_complex[i] * trace_real1[i]);
}
/* apply ifft */
FFT(-1,(unsigned long)(log(npadd)/log(2.0)),trace_real_conv,trace_complex_conv);
/* write output into data matrix */
for(j=1;j<=ns;j++){
temp_conv[j] = trace_real_conv[j];
}
free_vector(trace_real,1,npad);
free_vector(trace_complex,1,npad);
free_vector(trace_real_conv,1,npad);
free_vector(trace_complex_conv,1,npad);
free_vector(trace_real1,1,npad);
free_vector(trace_complex1,1,npad);
free_vector(window,1,(npad/2));
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2013 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* send L2 norm of each PE to PE 0
* ----------------------------------------------------------------------*/
#include "fd.h"
float exchange_L2(float L2,int sw, int bcast_l2){
extern int NP, MYID;
float l2, *L2s, *L2c;
int i;
const int tag=1;
MPI_Status status;
L2s=vector(1,2); /* vector for L2 exchange */
L2c=vector(0,NP+1); /* vector for merging of L2 norms */
L2s[1]=L2;
if(MYID==0){L2c[0]=L2;}
/* send L2 norm values to PE 0 */
for(i=1;i<=NP-1;i++){
/* if(MYID==i){
printf("Process %d sending %e to %d\n",i, L2s[1], 0);}*/
if(MYID==i){
MPI_Send(&L2s[1],2,MPI_FLOAT,0,tag,MPI_COMM_WORLD);}